{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bikeshare预测\n",
    "1) 请在 Capital Bikeshare （美国 Washington, D.C.的一个共享单车公司）提供的自行 车数据上进行回归分析。训练数据为 2011 年的数据，要求预测 2012 年每天的单车 共享数量\n",
    "\n",
    "原始数据集地址：http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset 1) 文件说明 day.csv: 按天计的单车共享次数（作业只需使用该文件） hour.csv: 按小时计的单车共享次数（无需理会） readme：数据说明文件\n",
    "\n",
    "2) 字段说明 Instant记录号 Dteday：日期 Season：季节（1=春天、2=夏天、3=秋天、4=冬天） yr：年份，(0: 2011, 1:2012) mnth：月份( 1 to 12) hr：小时 (0 to 23) （只在hour.csv有，作业忽略此字段） holiday：是否是节假日 weekday：星期中的哪天，取值为0～6 workingday：是否工作日 1=工作日 （是否为工作日，1为工作日，0为非周末或节假日 weathersit：天气（1：晴天，多云 ",
    "2：雾天，阴天 ",
    "3：小雪，小雨 ",
    "4：大雨，大雪，大雾） temp：气温摄氏度 atemp：体感温度 hum：湿度 windspeed：风速\n",
    "\n",
    "casual：非注册用户个数 registered：注册用户个数 cnt：给定日期（天）时间（每小时）总租车人数，响应变量y casual、registered和cnt三个特征均为要预测的y，作业里只需对cnt进行预测\n",
    "\n",
    "**********************************************************************************************************************************************\n",
    "季节，月份，星期，天气，是类别特征，用数字表示是不合适的，需要独热编码（One hot encoding）/哑编码\n",
    "是否是节假日，是否工作日，这样只有0，1的判别特征不需要标准化\n",
    "\n",
    "思路：\n",
    "1）读入数据：用pandas读入csv文件，根据yr（年份）标签的0，1值将数据分成训练集和测试集\n",
    "2）数据探索：大概看一下数据信息，分布，有没有空值等\n",
    "3）特征工程：机器学习还是需要人工处理特征，一般是对连续值特征标准化（也有MinMaxScaler等方法），类别特征做OneHot编码等，具体看特征工程ppt\n",
    "4）确定模型：本题是回归问题，题目要求用Lasso和岭回归\n",
    "5）模型训练："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "#模型，ElasticNetCV（弹性网络）这里没用\n",
    "from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV\n",
    "\n",
    "#模型评估\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "#可视化\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>dteday</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>casual</th>\n",
       "      <th>registered</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2011-01-01</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.344167</td>\n",
       "      <td>0.363625</td>\n",
       "      <td>0.805833</td>\n",
       "      <td>0.160446</td>\n",
       "      <td>331</td>\n",
       "      <td>654</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>2011-01-02</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.363478</td>\n",
       "      <td>0.353739</td>\n",
       "      <td>0.696087</td>\n",
       "      <td>0.248539</td>\n",
       "      <td>131</td>\n",
       "      <td>670</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>2011-01-03</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.196364</td>\n",
       "      <td>0.189405</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.248309</td>\n",
       "      <td>120</td>\n",
       "      <td>1229</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2011-01-04</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.212122</td>\n",
       "      <td>0.590435</td>\n",
       "      <td>0.160296</td>\n",
       "      <td>108</td>\n",
       "      <td>1454</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>2011-01-05</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.226957</td>\n",
       "      <td>0.229270</td>\n",
       "      <td>0.436957</td>\n",
       "      <td>0.186900</td>\n",
       "      <td>82</td>\n",
       "      <td>1518</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant      dteday  season  yr  mnth  holiday  weekday  workingday  \\\n",
       "0        1  2011-01-01       1   0     1        0        6           0   \n",
       "1        2  2011-01-02       1   0     1        0        0           0   \n",
       "2        3  2011-01-03       1   0     1        0        1           1   \n",
       "3        4  2011-01-04       1   0     1        0        2           1   \n",
       "4        5  2011-01-05       1   0     1        0        3           1   \n",
       "\n",
       "   weathersit      temp     atemp       hum  windspeed  casual  registered  \\\n",
       "0           2  0.344167  0.363625  0.805833   0.160446     331         654   \n",
       "1           2  0.363478  0.353739  0.696087   0.248539     131         670   \n",
       "2           1  0.196364  0.189405  0.437273   0.248309     120        1229   \n",
       "3           1  0.200000  0.212122  0.590435   0.160296     108        1454   \n",
       "4           1  0.226957  0.229270  0.436957   0.186900      82        1518   \n",
       "\n",
       "    cnt  \n",
       "0   985  \n",
       "1   801  \n",
       "2  1349  \n",
       "3  1562  \n",
       "4  1600  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#读入数据\n",
    "data = pd.read_csv(\"day.csv\")\n",
    "data.head()\n",
    "#print(\"train : \" + str(train.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 731 entries, 0 to 730\n",
      "Data columns (total 16 columns):\n",
      "instant       731 non-null int64\n",
      "dteday        731 non-null object\n",
      "season        731 non-null int64\n",
      "yr            731 non-null int64\n",
      "mnth          731 non-null int64\n",
      "holiday       731 non-null int64\n",
      "weekday       731 non-null int64\n",
      "workingday    731 non-null int64\n",
      "weathersit    731 non-null int64\n",
      "temp          731 non-null float64\n",
      "atemp         731 non-null float64\n",
      "hum           731 non-null float64\n",
      "windspeed     731 non-null float64\n",
      "casual        731 non-null int64\n",
      "registered    731 non-null int64\n",
      "cnt           731 non-null int64\n",
      "dtypes: float64(4), int64(11), object(1)\n",
      "memory usage: 91.5+ KB\n"
     ]
    }
   ],
   "source": [
    "#看有没有缺失值，各特征数据类型，行列数\n",
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据探索"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>casual</th>\n",
       "      <th>registered</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "      <td>731.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>366.000000</td>\n",
       "      <td>2.496580</td>\n",
       "      <td>0.500684</td>\n",
       "      <td>6.519836</td>\n",
       "      <td>0.028728</td>\n",
       "      <td>2.997264</td>\n",
       "      <td>0.683995</td>\n",
       "      <td>1.395349</td>\n",
       "      <td>0.495385</td>\n",
       "      <td>0.474354</td>\n",
       "      <td>0.627894</td>\n",
       "      <td>0.190486</td>\n",
       "      <td>848.176471</td>\n",
       "      <td>3656.172367</td>\n",
       "      <td>4504.348837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>211.165812</td>\n",
       "      <td>1.110807</td>\n",
       "      <td>0.500342</td>\n",
       "      <td>3.451913</td>\n",
       "      <td>0.167155</td>\n",
       "      <td>2.004787</td>\n",
       "      <td>0.465233</td>\n",
       "      <td>0.544894</td>\n",
       "      <td>0.183051</td>\n",
       "      <td>0.162961</td>\n",
       "      <td>0.142429</td>\n",
       "      <td>0.077498</td>\n",
       "      <td>686.622488</td>\n",
       "      <td>1560.256377</td>\n",
       "      <td>1937.211452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.059130</td>\n",
       "      <td>0.079070</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.022392</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>20.000000</td>\n",
       "      <td>22.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>183.500000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.337083</td>\n",
       "      <td>0.337842</td>\n",
       "      <td>0.520000</td>\n",
       "      <td>0.134950</td>\n",
       "      <td>315.500000</td>\n",
       "      <td>2497.000000</td>\n",
       "      <td>3152.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>366.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>7.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.498333</td>\n",
       "      <td>0.486733</td>\n",
       "      <td>0.626667</td>\n",
       "      <td>0.180975</td>\n",
       "      <td>713.000000</td>\n",
       "      <td>3662.000000</td>\n",
       "      <td>4548.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>548.500000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>10.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>0.655417</td>\n",
       "      <td>0.608602</td>\n",
       "      <td>0.730209</td>\n",
       "      <td>0.233214</td>\n",
       "      <td>1096.000000</td>\n",
       "      <td>4776.500000</td>\n",
       "      <td>5956.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>731.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>12.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>0.861667</td>\n",
       "      <td>0.840896</td>\n",
       "      <td>0.972500</td>\n",
       "      <td>0.507463</td>\n",
       "      <td>3410.000000</td>\n",
       "      <td>6946.000000</td>\n",
       "      <td>8714.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          instant      season          yr        mnth     holiday     weekday  \\\n",
       "count  731.000000  731.000000  731.000000  731.000000  731.000000  731.000000   \n",
       "mean   366.000000    2.496580    0.500684    6.519836    0.028728    2.997264   \n",
       "std    211.165812    1.110807    0.500342    3.451913    0.167155    2.004787   \n",
       "min      1.000000    1.000000    0.000000    1.000000    0.000000    0.000000   \n",
       "25%    183.500000    2.000000    0.000000    4.000000    0.000000    1.000000   \n",
       "50%    366.000000    3.000000    1.000000    7.000000    0.000000    3.000000   \n",
       "75%    548.500000    3.000000    1.000000   10.000000    0.000000    5.000000   \n",
       "max    731.000000    4.000000    1.000000   12.000000    1.000000    6.000000   \n",
       "\n",
       "       workingday  weathersit        temp       atemp         hum   windspeed  \\\n",
       "count  731.000000  731.000000  731.000000  731.000000  731.000000  731.000000   \n",
       "mean     0.683995    1.395349    0.495385    0.474354    0.627894    0.190486   \n",
       "std      0.465233    0.544894    0.183051    0.162961    0.142429    0.077498   \n",
       "min      0.000000    1.000000    0.059130    0.079070    0.000000    0.022392   \n",
       "25%      0.000000    1.000000    0.337083    0.337842    0.520000    0.134950   \n",
       "50%      1.000000    1.000000    0.498333    0.486733    0.626667    0.180975   \n",
       "75%      1.000000    2.000000    0.655417    0.608602    0.730209    0.233214   \n",
       "max      1.000000    3.000000    0.861667    0.840896    0.972500    0.507463   \n",
       "\n",
       "            casual   registered          cnt  \n",
       "count   731.000000   731.000000   731.000000  \n",
       "mean    848.176471  3656.172367  4504.348837  \n",
       "std     686.622488  1560.256377  1937.211452  \n",
       "min       2.000000    20.000000    22.000000  \n",
       "25%     315.500000  2497.000000  3152.000000  \n",
       "50%     713.000000  3662.000000  4548.000000  \n",
       "75%    1096.000000  4776.500000  5956.000000  \n",
       "max    3410.000000  6946.000000  8714.000000  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#对数据值型特征，用常用统计量观察其分布\n",
    "#count(总个数)， mean(均值)，std(标准差)，分位数，最大最小值\n",
    "data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "season属性的不同取值和出现的次数\n",
      "3    188\n",
      "2    184\n",
      "1    181\n",
      "4    178\n",
      "Name: season, dtype: int64\n",
      "\n",
      "mnth属性的不同取值和出现的次数\n",
      "12    62\n",
      "10    62\n",
      "8     62\n",
      "7     62\n",
      "5     62\n",
      "3     62\n",
      "1     62\n",
      "11    60\n",
      "9     60\n",
      "6     60\n",
      "4     60\n",
      "2     57\n",
      "Name: mnth, dtype: int64\n",
      "\n",
      "weathersit属性的不同取值和出现的次数\n",
      "1    463\n",
      "2    247\n",
      "3     21\n",
      "Name: weathersit, dtype: int64\n",
      "\n",
      "weekday属性的不同取值和出现的次数\n",
      "6    105\n",
      "1    105\n",
      "0    105\n",
      "5    104\n",
      "4    104\n",
      "3    104\n",
      "2    104\n",
      "Name: weekday, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "#对类别型特征，观察其取值范围及直方图\n",
    "#   value_counts()   对数据统计，返回 出现的值-频次 的对儿\n",
    "#   astype()    强制类型转换，这里变成了object类型\n",
    "categorical_features = ['season','mnth','weathersit','weekday']\n",
    "for col in categorical_features:\n",
    "    print ('\\n%s属性的不同取值和出现的次数'%col)\n",
    "    print (data[col].value_counts())\n",
    "    data[col] = data[col].astype('object')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 731 entries, 0 to 730\n",
      "Data columns (total 16 columns):\n",
      "instant       731 non-null int64\n",
      "dteday        731 non-null object\n",
      "season        731 non-null object\n",
      "yr            731 non-null int64\n",
      "mnth          731 non-null object\n",
      "holiday       731 non-null int64\n",
      "weekday       731 non-null object\n",
      "workingday    731 non-null int64\n",
      "weathersit    731 non-null object\n",
      "temp          731 non-null float64\n",
      "atemp         731 non-null float64\n",
      "hum           731 non-null float64\n",
      "windspeed     731 non-null float64\n",
      "casual        731 non-null int64\n",
      "registered    731 non-null int64\n",
      "cnt           731 non-null int64\n",
      "dtypes: float64(4), int64(7), object(5)\n",
      "memory usage: 91.5+ KB\n"
     ]
    }
   ],
   "source": [
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 特征工程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>season_1</th>\n",
       "      <th>season_2</th>\n",
       "      <th>season_3</th>\n",
       "      <th>season_4</th>\n",
       "      <th>mnth_1</th>\n",
       "      <th>mnth_2</th>\n",
       "      <th>mnth_3</th>\n",
       "      <th>mnth_4</th>\n",
       "      <th>mnth_5</th>\n",
       "      <th>mnth_6</th>\n",
       "      <th>...</th>\n",
       "      <th>weathersit_1</th>\n",
       "      <th>weathersit_2</th>\n",
       "      <th>weathersit_3</th>\n",
       "      <th>weekday_0</th>\n",
       "      <th>weekday_1</th>\n",
       "      <th>weekday_2</th>\n",
       "      <th>weekday_3</th>\n",
       "      <th>weekday_4</th>\n",
       "      <th>weekday_5</th>\n",
       "      <th>weekday_6</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  mnth_4  \\\n",
       "0         1         0         0         0       1       0       0       0   \n",
       "1         1         0         0         0       1       0       0       0   \n",
       "2         1         0         0         0       1       0       0       0   \n",
       "3         1         0         0         0       1       0       0       0   \n",
       "4         1         0         0         0       1       0       0       0   \n",
       "\n",
       "   mnth_5  mnth_6    ...      weathersit_1  weathersit_2  weathersit_3  \\\n",
       "0       0       0    ...                 0             1             0   \n",
       "1       0       0    ...                 0             1             0   \n",
       "2       0       0    ...                 1             0             0   \n",
       "3       0       0    ...                 1             0             0   \n",
       "4       0       0    ...                 1             0             0   \n",
       "\n",
       "   weekday_0  weekday_1  weekday_2  weekday_3  weekday_4  weekday_5  weekday_6  \n",
       "0          0          0          0          0          0          0          1  \n",
       "1          1          0          0          0          0          0          0  \n",
       "2          0          1          0          0          0          0          0  \n",
       "3          0          0          1          0          0          0          0  \n",
       "4          0          0          0          1          0          0          0  \n",
       "\n",
       "[5 rows x 26 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#类别特征OneHot编码\n",
    "#   pd.get_dummies()   OneHot编码函数\n",
    "categorical_features = ['season','mnth','weathersit','weekday']\n",
    "cat_features_block = data[categorical_features]\n",
    "cat_features_block = pd.get_dummies(cat_features_block)\n",
    "cat_features_block.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.826662</td>\n",
       "      <td>-0.679946</td>\n",
       "      <td>1.250171</td>\n",
       "      <td>-0.387892</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.721095</td>\n",
       "      <td>-0.740652</td>\n",
       "      <td>0.479113</td>\n",
       "      <td>0.749602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-1.634657</td>\n",
       "      <td>-1.749767</td>\n",
       "      <td>-1.339274</td>\n",
       "      <td>0.746632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-1.614780</td>\n",
       "      <td>-1.610270</td>\n",
       "      <td>-0.263182</td>\n",
       "      <td>-0.389829</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.467414</td>\n",
       "      <td>-1.504971</td>\n",
       "      <td>-1.341494</td>\n",
       "      <td>-0.046307</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       temp     atemp       hum  windspeed\n",
       "0 -0.826662 -0.679946  1.250171  -0.387892\n",
       "1 -0.721095 -0.740652  0.479113   0.749602\n",
       "2 -1.634657 -1.749767 -1.339274   0.746632\n",
       "3 -1.614780 -1.610270 -0.263182  -0.389829\n",
       "4 -1.467414 -1.504971 -1.341494  -0.046307"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#   fit_transform()函数返回值是numpy.ndarray类型\n",
    "#   后面三种线性回归模型里的训练函数fit()，numpy.ndarray和pandas的DataFram两种类型都能接受\n",
    "#数值型特征预处理（标准化等方法）\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "#标准化器\n",
    "ss_X = StandardScaler()\n",
    "\n",
    "numerical_features = ['temp','atemp','hum','windspeed']\n",
    "temp = ss_X.fit_transform(data[numerical_features])\n",
    "\n",
    "num_featurea_block = pd.DataFrame(data=temp, columns=numerical_features, index=data.index)\n",
    "num_featurea_block.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>season_1</th>\n",
       "      <th>season_2</th>\n",
       "      <th>season_3</th>\n",
       "      <th>season_4</th>\n",
       "      <th>mnth_1</th>\n",
       "      <th>mnth_2</th>\n",
       "      <th>mnth_3</th>\n",
       "      <th>mnth_4</th>\n",
       "      <th>mnth_5</th>\n",
       "      <th>mnth_6</th>\n",
       "      <th>...</th>\n",
       "      <th>weekday_3</th>\n",
       "      <th>weekday_4</th>\n",
       "      <th>weekday_5</th>\n",
       "      <th>weekday_6</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>holiday</th>\n",
       "      <th>workingday</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.826662</td>\n",
       "      <td>-0.679946</td>\n",
       "      <td>1.250171</td>\n",
       "      <td>-0.387892</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.721095</td>\n",
       "      <td>-0.740652</td>\n",
       "      <td>0.479113</td>\n",
       "      <td>0.749602</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.634657</td>\n",
       "      <td>-1.749767</td>\n",
       "      <td>-1.339274</td>\n",
       "      <td>0.746632</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.614780</td>\n",
       "      <td>-1.610270</td>\n",
       "      <td>-0.263182</td>\n",
       "      <td>-0.389829</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.467414</td>\n",
       "      <td>-1.504971</td>\n",
       "      <td>-1.341494</td>\n",
       "      <td>-0.046307</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 32 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  mnth_4  \\\n",
       "0         1         0         0         0       1       0       0       0   \n",
       "1         1         0         0         0       1       0       0       0   \n",
       "2         1         0         0         0       1       0       0       0   \n",
       "3         1         0         0         0       1       0       0       0   \n",
       "4         1         0         0         0       1       0       0       0   \n",
       "\n",
       "   mnth_5  mnth_6     ...      weekday_3  weekday_4  weekday_5  weekday_6  \\\n",
       "0       0       0     ...              0          0          0          1   \n",
       "1       0       0     ...              0          0          0          0   \n",
       "2       0       0     ...              0          0          0          0   \n",
       "3       0       0     ...              0          0          0          0   \n",
       "4       0       0     ...              1          0          0          0   \n",
       "\n",
       "       temp     atemp       hum  windspeed  holiday  workingday  \n",
       "0 -0.826662 -0.679946  1.250171  -0.387892        0           0  \n",
       "1 -0.721095 -0.740652  0.479113   0.749602        0           0  \n",
       "2 -1.634657 -1.749767 -1.339274   0.746632        0           1  \n",
       "3 -1.614780 -1.610270 -0.263182  -0.389829        0           1  \n",
       "4 -1.467414 -1.504971 -1.341494  -0.046307        0           1  \n",
       "\n",
       "[5 rows x 32 columns]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#   pd.concat   合并行/列的函数，这里合并列令axis=1，ignore_index=False不忽略index数据根据index对齐\n",
    "#   合并出来的列顺序跟pd.concat里的列名顺序一样\n",
    "x_data = pd.concat([cat_features_block, num_featurea_block, data['holiday'], data['workingday']], axis=1, ignore_index=False)\n",
    "x_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season_1</th>\n",
       "      <th>season_2</th>\n",
       "      <th>season_3</th>\n",
       "      <th>season_4</th>\n",
       "      <th>mnth_1</th>\n",
       "      <th>mnth_2</th>\n",
       "      <th>mnth_3</th>\n",
       "      <th>mnth_4</th>\n",
       "      <th>mnth_5</th>\n",
       "      <th>...</th>\n",
       "      <th>weekday_5</th>\n",
       "      <th>weekday_6</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>holiday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>yr</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.826662</td>\n",
       "      <td>-0.679946</td>\n",
       "      <td>1.250171</td>\n",
       "      <td>-0.387892</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.721095</td>\n",
       "      <td>-0.740652</td>\n",
       "      <td>0.479113</td>\n",
       "      <td>0.749602</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.634657</td>\n",
       "      <td>-1.749767</td>\n",
       "      <td>-1.339274</td>\n",
       "      <td>0.746632</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.614780</td>\n",
       "      <td>-1.610270</td>\n",
       "      <td>-0.263182</td>\n",
       "      <td>-0.389829</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.467414</td>\n",
       "      <td>-1.504971</td>\n",
       "      <td>-1.341494</td>\n",
       "      <td>-0.046307</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 35 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant  season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  \\\n",
       "0        1         1         0         0         0       1       0       0   \n",
       "1        2         1         0         0         0       1       0       0   \n",
       "2        3         1         0         0         0       1       0       0   \n",
       "3        4         1         0         0         0       1       0       0   \n",
       "4        5         1         0         0         0       1       0       0   \n",
       "\n",
       "   mnth_4  mnth_5  ...   weekday_5  weekday_6      temp     atemp       hum  \\\n",
       "0       0       0  ...           0          1 -0.826662 -0.679946  1.250171   \n",
       "1       0       0  ...           0          0 -0.721095 -0.740652  0.479113   \n",
       "2       0       0  ...           0          0 -1.634657 -1.749767 -1.339274   \n",
       "3       0       0  ...           0          0 -1.614780 -1.610270 -0.263182   \n",
       "4       0       0  ...           0          0 -1.467414 -1.504971 -1.341494   \n",
       "\n",
       "   windspeed  holiday  workingday  yr   cnt  \n",
       "0  -0.387892        0           0   0   985  \n",
       "1   0.749602        0           0   0   801  \n",
       "2   0.746632        0           1   0  1349  \n",
       "3  -0.389829        0           1   0  1562  \n",
       "4  -0.046307        0           1   0  1600  \n",
       "\n",
       "[5 rows x 35 columns]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#把剩下的列补上生产完整的新表格，是个好习惯\n",
    "FE_data = pd.concat([data['instant'], x_data, data['yr'], data['cnt']], axis=1)\n",
    "FE_data.to_csv('FE_data.csv', index=False)\n",
    "FE_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 731 entries, 0 to 730\n",
      "Data columns (total 35 columns):\n",
      "instant         731 non-null int64\n",
      "season_1        731 non-null uint8\n",
      "season_2        731 non-null uint8\n",
      "season_3        731 non-null uint8\n",
      "season_4        731 non-null uint8\n",
      "mnth_1          731 non-null uint8\n",
      "mnth_2          731 non-null uint8\n",
      "mnth_3          731 non-null uint8\n",
      "mnth_4          731 non-null uint8\n",
      "mnth_5          731 non-null uint8\n",
      "mnth_6          731 non-null uint8\n",
      "mnth_7          731 non-null uint8\n",
      "mnth_8          731 non-null uint8\n",
      "mnth_9          731 non-null uint8\n",
      "mnth_10         731 non-null uint8\n",
      "mnth_11         731 non-null uint8\n",
      "mnth_12         731 non-null uint8\n",
      "weathersit_1    731 non-null uint8\n",
      "weathersit_2    731 non-null uint8\n",
      "weathersit_3    731 non-null uint8\n",
      "weekday_0       731 non-null uint8\n",
      "weekday_1       731 non-null uint8\n",
      "weekday_2       731 non-null uint8\n",
      "weekday_3       731 non-null uint8\n",
      "weekday_4       731 non-null uint8\n",
      "weekday_5       731 non-null uint8\n",
      "weekday_6       731 non-null uint8\n",
      "temp            731 non-null float64\n",
      "atemp           731 non-null float64\n",
      "hum             731 non-null float64\n",
      "windspeed       731 non-null float64\n",
      "holiday         731 non-null int64\n",
      "workingday      731 non-null int64\n",
      "yr              731 non-null int64\n",
      "cnt             731 non-null int64\n",
      "dtypes: float64(4), int64(5), uint8(26)\n",
      "memory usage: 70.0 KB\n"
     ]
    }
   ],
   "source": [
    "FE_data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "至此数据预处理完成"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bike Sharing数据集上的回归分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(731, 35)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#读入数据\n",
    "data = pd.read_csv(\"FE_data.csv\")\n",
    "data.head()\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "reset_index(drop=True)  重置index，最后保持测试集原来的index比较好，这里不用重置index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train :(365, 33)\n",
      "test :(366, 33)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:4: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  after removing the cwd from sys.path.\n",
      "E:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:9: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  if __name__ == '__main__':\n"
     ]
    }
   ],
   "source": [
    "#2011年做训练集，2012年做测试集\n",
    "#   reset_index(drop=True)  重置index\n",
    "train = data[data.yr == 0]\n",
    "train.drop(['instant', 'yr'], axis=1, inplace=True)\n",
    "print(\"train :\" + str(train.shape))\n",
    "\n",
    "test = data[data.yr == 1]\n",
    "testID = test['instant']\n",
    "test.drop(['instant', 'yr'], axis=1, inplace=True)\n",
    "#test = test.reset_index(drop=True)\n",
    "print(\"test :\" + str(test.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备训练数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = train[\"cnt\"]\n",
    "x_train = train.drop(['cnt'], axis=1)\n",
    "\n",
    "y_test = test[\"cnt\"]\n",
    "x_test = test.drop(['cnt'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean y diff : 1.5914175510574313\n"
     ]
    }
   ],
   "source": [
    "#对y标准化\n",
    "#可以在特征工程哪里一起做了\n",
    "#====对y标准化不是必须，但对其进行标准化可以使得不同问题w的取值范围相对相同===#\n",
    "#这些参数需要保留，对测试集预测完后还需要对其进行反变换\n",
    "\n",
    "y_mean = y_train.mean()\n",
    "y_std = y_train.std()\n",
    "y_train = (y_train-y_mean)/y_std\n",
    "\n",
    "y_test = (y_test-y_mean)/y_std\n",
    "\n",
    "mean_y_test = y_test.mean()\n",
    "mean_diff = mean_y_test\n",
    "\n",
    "print(\"mean y diff :\", mean_diff)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 Linear Regression with Ridge regularization (L2 penalty)  岭回归L2正则"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于本数据集2011和2012年差距较大加入了均值补偿\n",
    "y_test_pred += mean_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 1.0\n",
      "cv of rmse : 0.4354843278524383\n",
      "RMSE on Training set : 0.3991636934921292\n",
      "RMSE on Test set : 0.7438808968273131\n",
      "r2_score on Training set : 0.84023062147299\n",
      "r2_score on Test set : 0.6703075644943807\n"
     ]
    }
   ],
   "source": [
    "#岭回归／L2正则\n",
    "#   scoring(评价指标，缺省值默认MSE)， cv（交叉验证，默认留1交叉验证），gcv_mode（留1交叉验证的模式）\n",
    "#   store_cv_values（是否保存cv的值）\n",
    "#class sklearn.linear_model.RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, \n",
    "#                                  normalize=False, scoring=None, cv=None, gcv_mode=None, \n",
    "#                                  store_cv_values=False)\n",
    "\n",
    "#生成学习器实例，这里自定义了alphas\n",
    "alphas = [0.01, 0.1, 1, 10, 100]\n",
    "ridge = RidgeCV(alphas=alphas, store_cv_values=True)\n",
    "\n",
    "#训练\n",
    "#RidgeCV采用的是广义交叉验证（Generalized Cross-Validation），留一交叉验证（N-折交叉验证）的一种有效实现方式\n",
    "#RidgeCV封装了交叉验证，不用定义校验集了\n",
    "ridge.fit(x_train, y_train)\n",
    "\n",
    "#通过交叉验证得到的最佳超参数alpha\n",
    "alpha = ridge.alpha_\n",
    "print(\"Best alpha :\", alpha)\n",
    "\n",
    "# 交叉验证估计的测试误差\n",
    "mse_cv = np.mean(ridge.cv_values_, axis = 0)\n",
    "rmse_cv = np.sqrt(mse_cv)\n",
    "print(\"cv of rmse :\", min(rmse_cv))\n",
    "\n",
    "#训练上测试，训练误差，实际任务中这一步不需要\n",
    "y_train_pred = ridge.predict(x_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "\n",
    "y_test_pred = ridge.predict(x_test)\n",
    "y_test_pred += mean_diff\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\", rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\", r2_score_train)\n",
    "print(\"r2_score on Test set :\", r2_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ridge picked 32 features and eliminated the other 0 features\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEICAYAAAD7pTujAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcXFWd/vHPIzsEiUpkQIRogImyBQkgP0ABERzcRxABlTBuiA7jOCAMOpgRHVHccEEERBBQEQSGTQFhIDGGsGYBQRCCQwTZFMQRIpjn98c5TSpNL1XdnarqzvN+vfqVqnvPPfd7qjt9+px77vfKNhEREd3keZ0OICIiord0ThER0XXSOUVERNdJ5xQREV0nnVNERHSddE4REdF10jnFqCbp7yXdIukJSYdJWkPSxZIel3SupAMlXdFEPUdLOrUdMQ8Qw0aS/ixppRGqb7qks0airibONeDnLOkaSe9vRywjoZV4JVnSJss7phVNOqdoC0kHSLqx/vJ9QNJPJe08AlV/ArjG9tq2vw7sA6wHvMj2vrbPtr3nYJXY/i/bw/7lKWli/WW1cqvH2v5f2+Ns/20I591V0qJWj2uh/tMl/bV+//4g6UpJk3v2N/s5L4e4ptfP+7Be2z9Wt09vd0wxMtI5xXIn6ePA14D/onQcGwEnAm8dgeo3Bm7r9f5O28+MQN2xrC/aHge8BPgd8N0Ox9PjTuCgXtveW7fHKJXOKZYrSesAnwE+Yvt82/9n+2nbF9s+opZZTdLXJN1fv74mabWGOt4kaa6kxyT9UtJWdfvVwG7AN+tf9D8EjgH2q+/fJ2mapF801LV5/av/D5IelHR03b7MFJikV9dzPSZpnqRdG/ZdI+lYSbPqdOIVktatu2fUfx+rMewoaRNJ19apxkckndPPZ7XMqGuQ8zQetxbwU2CDes4/S9qg7l5V0vfr8bdJmtpw3AaSfiLpYUkLe48++mP7SeDHwJSGunp/zq+XdEdt8zcBNexbSdKX62exUNJHe7V7HUnfrSPs30n67CBTnTcAa0ravB6/ObBG3d74OX1A0m/q9/6ihs9owHjr/n+SdLukP0q6XNLGzXxWMXTpnGJ52xFYHbhggDKfBF5N+WW3NbA98CkASa8CTgM+BLwI+A5wkaTVbO8OzAQ+WqfD9qeMzs6p75f5y17S2sDPgZ8BGwCbAFf1DkbSS4BLgc8CLwQOB34iaUJDsQOAg4EXA6vWMgCvqf+OrzHMBo4FrgBeAGwIfGOAz6K3/s7zLNv/B/wDcH895zjb99fdbwF+BIwHLgK+Wdv4POBiYB5lJPQ64GOS9hosoNoZ7g/8pp/96wI/oXwP1wXuBnZqKPKBGu8U4FXA23pVcQbwDOX7sw2wJzDYlOuZlNESlFHU93vFtDvweeCdwPrAbymfy6DxSnobcDTwj8AEys/cDweJJ4YpnVMsby8CHhlkmu1A4DO2H7L9MPCfwHvqvg8A37E9x/bfbJ8BLKZ0Zq16E/B721+2/ZTtJ2zP6aPcu4HLbF9me4ntK4Ebgb0bynzP9p19jSL68DRlunGDet5fDFC2t1bO05df1Hb8jfILfOu6fTtggu3P2P6r7XuAU4B3DVDX4ZIeA54Admbp96i3vYFf2T7P9tOUKd3fN+x/J3CC7UW2/wgc17ND0nqUjutjdZT9EPDVQeICOAvYX9IqtWzvhSAHAqfZvtn2YuDfgR0lTWwi3g8Bn7d9e/05/i9gSkZPy1c6p1jeHgXW1cALBDag/CXb47d1G5Rf6v9Wp9ceq78cX9qwvxUvpfxVPJiNgX17nXNnyl/cPRp/ef0FGDdAfZ+gTBNdX6fW/qmFmFs5TzPHr16/FxtTpgEb23g05Zpgf75kezwwEXgS+Pt+ym0A3NfzxiW79H397e/1emNgFeCBhri+Qxk59sv2/1JGcv8F3GX7vl5FlvkZs/1nys/mS5qId2PghIZ4/kD5fr5koJhieFpeURTRotnAU5Spm/P6KXM/yy5s2Khug/JL4nO2PzcCsdxHmY5qptyZtj8whHM8J82/7d9TRoCorFD8uaQZtvucFhuiVh8vcB+w0PamLZ/I/l9J/wKcIemSOqpr9ADlDwEAJKnxfd2/YcP7xn33UUbG6w5hUcv3KVPAB/exr+dnrCemtSij+t81EW/Pz+DZLcYTw5CRUyxXth+nLFL4lqS3SVpT0iqS/kHSF2uxHwKfkjShzv8fw9JpmVOAQyTtoGItSW+s149adQnwdyrLjFeTtLakHfoodxbwZkl71Yv3q6ss1d6wj7K9PQwsAV7es0HSvg3H/pHSkbS8XHwQDwIvUlmA0ozrgT9JOlLl3rCVJG0habtmDq5TnfcDH+xj96XA5pL+sY7SDgP+rmH/j4F/kfQSSeOBIxvqfYByfe7Lkp4v6XmSJkl6bRNhnUO5PvXjPvb9ADhY0hSVxTb/BcyxfW8T8Z4E/HvDgot1JO3bRDwxDOmcYrmz/RXg45QLzg9T/hL9KHBhLfJZyjWd+cAC4Oa6Dds3UkYd36T8Yv8NMG2IcTwBvB54M2W66y7Kar/e5e6jLHM/uiHeI2ji/4vtvwCfA2bVaaBXU67vzJH0Z8qihH+xvXAobRjgvHdQOvl76nkHnPas16DeTLmGtRB4BDgVaLZzAzge+IQaVlbWuh8B9qVcS3oU2BSY1VDkFEoHNB+4BbiMsgCip8N+L2Xxx68o3/PzWHZKtb82PWn7532M5LB9FfAflIUPDwCTqNexBovX9gXAF4AfSfoTcCvlulgsR8rDBiOikyT9A3CS7SwwiGdl5BQRbVWnEfeWtHJdtv9pBr7VIFZAGTlFRFtJWhO4FphMWfV3KWWq808dDSy6SjqniIjoOpnWi4iIrpP7nIZo3XXX9cSJEzsdRkTEqHLTTTc9YnvCYOXSOQ3RxIkTufHGGzsdRkTEqCLpt4OXyrReRER0oYycIgKAq66e1OkQYpR43e7NpKgcnjEzcpI0XtKhDe93lXRJC8efrvJsmbn1q9XszxERMULGTOdEeV7NoYOWGtgRtqfUr7kjEVRERLSuqzonlSeB3iHpVEm3Sjpb0h4qTwK9S9L2Kk8sPU3lKaH3aOnTO48DJtVRz/F12zhJ59U6z67ZhocT3wcl3SjpxocffnhYbY2IiP51VedUbQKcAGxFuYP8AMqzdA6nJOKkbt+L8sTUT9cHjB0F3F1HPUfUctsAHwNeSckS3fg0zr58TtJ8SV/tncwSwPbJtqfanjphwqArISMiYoi6sXNaaHuB7SWU5/tcVR/+tYDykDOAS20vrtmEH6L/B6RdX5+2uQSY23B8X/6d0ultR3k095EDlI2IiOWoG1frLW54vaTh/RKWxttY5m/0345my/U8RwZgsaTvUUZqESuMdqzAimhWN46chuoJYCgPoANA0vr1X1Ge2nrrCMUVEREt6saR05DYfrQunLgV+Ckl03ErzpY0ARBlCvCQkY4xIiKak6zkQzR16lQnfVFERGsk3WR76mDlxtK0XkREjBFtn9aTNA24wvb99f29wNS68m4kz3MZZRk6wAG2T6zbLwBe1qv4kcC/Aq8GfmH7TSMZS8RoMH369E6HEC0Y69+vTlxzmkZZbHD/cCuStLLtZ/raZ3vvWmYiJXPEiXX72/up6xlgTeBDw40rIiKGZ9BpPUmf6MnCUG9Ovbq+fp2ksyTtKWm2pJslnStpXN1/jKQbaqaHk1XsA0ylLD6YK2mNepp/rscvkDS5Hr9WzQRxg6RbJL21bp9Wz3MxcIWk9SXNqPXdKmmXWu5eSevSd+aI57B9FWXFX0REdFgz15xmALvU11MpKYFWoWRtWAB8CtjD9quAG4GP17LftL2d7S2ANYA32T6vljmwZnJ4spZ9pB7/bZbeX/RJ4Grb2wG7AcdLWqvu2xE4yPbulKm7y21PAbamrLRr1FfmiCFJ+qKIiPZopnO6CdhW0tqUm1pnUzqpXYAnKamBZkmaCxwEbFyP203SHEkLgN2BzQc4x/kN55pYX+8JHFXrvQZYHdio7rvS9h/q6xuAgyVNB7a0vdxGP0lfFBHRHoNec7L9dF20cDDwS2A+ZSQzCVhI6Sj2bzxG0uqUazxTbd9XO47VBzhNTyaHxiwOAt5h+9e96t4B+L+G+GZIeg3wRuBMScfb/v5g7YqIZY31C+wxujS7lHwGZbptBjCTcoPqXOA6YCdJmwBIWlPSZiztiB6p16D2aair2UwOl1OuRanWvU1fhSRtDDxk+xTgu8CrehUZVuaIiIhov2Y7p5nA+sBs2w8CTwEzbT9MWX33Q0nzKZ3VZNuPAadQrkldSJl663E6cFKvBRF9ORZYBZhfsz4c20+5XYG5km4B3kHJaP4s249Sph1vHWhBhKSZwLnA6yQtkrTXALFFRMRylAwRQ5QMERERrUuGiIiIGLVGbeJXSeNpyPzQ5DFbAmf22rzY9g4jGlxERAzLqO2cgPE0ZH5ohu0FwJTlFlHEKLboqJmdDmGFseFxuwxeaAU3mqf1lsn8IOmImk1ivqT/hJK6SNIdkk6tCyLOlrRHfbTGXZK2r+WmSzpT0tV1+wc62rKIiBXcaO6cns38AFwJbApsTxkZbVvvfQLYhLKCbyvKY9gPoGS3OBw4uqG+rSj3Su0IHCNpg3Y0IiIinms0d06N9qxftwA3UzqhTeu+hbYX2F4C3AZc5bJEcQFLs1EA/LftJ2t29P+hdHTLSPqiiIj2GM3XnBoJ+Lzt7yyzsWQkX9ywaUnD+yUs2/7ea+qfs8be9snAyVCWkg8r4oiI6Ndo7pwaMz9cDhwr6Wzbf5b0EuDpFut7q6TPA2tRbuw9asQijRgFcpE+usmo7ZxsP1oXNtwK/BT4ATC7Zjv6M/BuSq6+Zl0PXEpJLntsz8MQIyKi/UZt5wRg+4Bem07oo9gWDeWnNby+t3EfcKftD45kfBERMTRjZUFERESMIaN65DRSbE/vdAwREbFURk4REdF1MnLqRdJFwMvr4+UjVhhf3u9NnQ5hzPu3cy7pdAijRkZODST9I2WlX0REdFDbOidJa0m6VNK8muduP0nbSrpW0k2SLpe0fi37gZonb56kn0has27ftx47T9KMum11Sd+TtEDSLZJ2q9unSTpf0s9qvrwvDhLfOODjwGcHKJMMERERbdDOkdMbgPttb12nzH4GfAPYx/a2wGnA52rZ821vZ3tr4HbgfXX7McBedftb6raPANjeEtgfOENSz2PipwD7AVsC+0l66QDxHQt8GfhLfwVsn2x7qu2pEyZMaKXtERHRgnZ2TguAPSR9QdIuwEsp9xldKWku8Clgw1p2C0kzJS0ADgQ2r9tnAafXrOEr1W07U5/RZPsO4LfAZnXfVbYft/0U8Ctg474CkzQF2MT2BSPX3IiIGKq2LYiwfaekbYG9gc9TMonfZnvHPoqfDrzN9jxJ0yjphLB9iKQdKNnD59ZORQOctjGv3t/ov707UjKZ31vLvFjSNbZ3ba51EaNfLtZHN2nnNacNgL/YPgv4ErADMEHSjnX/KpJ6RkhrAw9IWoUycuqpY5LtObaPAR6hjL5m9JSRtBkl/dCvW4nN9rdtb2B7ImUkdmc6poiIzmnnUvItgeMlLaEkZf0w8AzwdUnr1Fi+RnmsxX8AcyhTdAtYmuD1eEmbUkZLVwHzgDuAk+oU4DPANNuLa469iIgYhVQebRStmjp1qm+88cZOhxERMapIusn21MHK5T6niIjoOitchghJc4DVem1+j+0FnYgnIiKea8x0TpLGAwfYPrG+3xU43PYyOVls79DP8aLcgLsvZWXft21/fbkGHdFFvnXI1Z0OYcz6yEm7dzqEUWcsTeuNBw4dxvHTKKv/Jtt+BfCjkQgqIiJa11Wdk6SJku6QdGpNU3S2pD3qE2/vkrS9pOmSTpN0jaR7JB1WDz8OmCRprqTj67Zxks6rdZ6tgZfwfRj4jO0lALYf6iO+pC+KiGiDruqcqk0oT7TdCpgMHEC59+hw4OhaZjKwF7A98Ol6P9RRwN22p9g+opbbBvgY8Erg5cBOA5x3EiXF0Y2SflqXrC8j6YsiItqjGzunhbYX1BHMbZQURKbc7zSxlrnU9mLbjwAPAev1U9f1thfVuuY2HN+X1YCn6hLHUyi5/iIiogO6cUFEY8qhJQ3vl7A03mbTEjVbDmAR8JP6+gLge80EGzFW5KJ9dJNuHDkN1RMszSQxFBcCPf87XwvcOeyIIiJiSLpx5DQkth+tCyduBX4KXNpiFccBZ0v6V8oDB98/0jFGRERzkr5oiJK+KCKidUlfFBERo9aYmdZrlqQLgJf12nyk7cs7EU9Et7h98is6HUJHveKO2zsdQjQYMyMnSeMlHdrwfldJz3l6mu2313uhGr8ul7S7pJvrzb9nSFrhOu6IiG4xZjonhpG+SNLzgDOAd9negvIcqYNGMLaIiGhBV3VOHUxf9CJgse2e5eNXAu/oI76kL4qIaIOu6pyqTqQvegRYRVLPCpJ9KElgl5H0RRER7dGNnVPb0xfV+t8FfFXS9ZQbep8ZqQZFRERruvGif0fSF9meDewCIGlPYLPmQ44Y/bJaLbpJN46chmpY6Yskvbj+uxpwJHDSCMUVEREtGjOdk+1HgVl1IcXxgx7wXEdIuh2YD1xsO48FjYjokKQvGqKkL4qIaF3SF0VExKjVjQsilqukL4ro25ZnbNnpENpiwUELOh1CNKHtIydJ0yRt0PD+XknrLofzXFZTGi2T1qiv9EXAg5JmS7pN0nxJ+410PBER0bxOTOtNAzYYrFAzBsp/Z3tv24/RXFqjvwDvtb058Abga5LGj0SMERHRukE7J0mf6EkRJOmrkq6ur18n6SxJe9ZRx82SzpU0ru4/RtINdfXcySr2AaZSHuo3V9Ia9TT/XI9fIGlyPX6tmqboBkm3SHpr3T6tnudi4ApJ60uaUeu7VVLPvUo9I7K+0hotw/adtu+qr++n3Nj7nBQQSV8UEdEezYycZlBvTqV0LONquqCdKVkbPgXsYftVwI3Ax2vZb9reriZSXQN4k+3zapkD65Tak7XsI/X4b1PSFAF8Erja9nbAbsDxktaq+3YEDrK9OyW90eV1em5rSiaIRn2lNeqXpO2BVYG7e+9L+qKIiPZoZkHETcC2ktamZFy4mdJJ7QJcRMlbN6vmVF0VmF2P203SJ4A1gRdSUhFd3M85zm841z/W13sCb5HU01mtDmxUX19p+w/19Q3AabXDvNB2786paZLWB86kdHxLhlpPREQMz6Cdk+2nJd0LHAz8knKT6m7AJGAhpaPYv/EYSasDJwJTbd8naTqlc+lPT5qhxhRDAt5h+9e96t4B+L+G+GZIeg3wRuBMScfb/v5g7epN0vOBS4FP2b6u1eMjRrusYotu0uyCiBmU6bYZwEzgEMr02XXATpI2AZC0pqTNWNoRPVKvQe3TUFezaYYup1yLUq17m74KSdoYeMj2KcB3gVf1KjLo+SStClwAfN/2uU3EFhERy1GzndNMYH1gtu0HgaeAmbYfpqy++6Gk+ZTOanJdJXcK5ZrUhZSptx6nAyf1WhDRl2OBVYD5km6t7/uyKzBX0i2UZzCd0LizybRG7wReA0yrcc2VNGWA2CIiYjlK+qIhSvqiiIjWJX1RRESMWitU+iJJW1JW4zVabHuHTsQT0VWmr9PpCIZm+uOdjiCWg1E1cuqdikjSrpIuaaGK1wLjKPdD7VHvfdqh3iD8dUm/qemLei+qiIiINhpVnRPNpSIayCxgD+C3vbb/A7Bp/fog5WbgiIjokE4kfp0o6Q5Jp9YVdGdL2kPSLEl3Sdpe0vSauugaSff0pE+i71RE4ySdV+s8u2fpeV9s32L73j52vZWyjNz1Hqfx9Ybc3rEnfVFERBt0auS0CWXJ91bAZEoKop0p91IdXctMBvYCtgc+XTNA9JWKaBvgY5RMFS8HdhpCPC8B7mt4v6huW0bSF0VEtEenOqeFthfUFEG3AVe5rGlfAEysZS61vdj2I5RErOv1U9f1thfVuuY2HN+KvkZbWWMfEdEhnVqtt7jh9ZKG90tYGlNjmca0RgPVNVC5gSwCXtrwfkPg/iHUEzF6ZdVbdJHRtiCi2dRHrboIeG9dtfdq4HHbDyyH80RERBNGVefUZCqifkk6TNIiyshovqRT667LgHuA31DSLg1nRWBERAxT0hcNUdIXRUS0LumLIiJi1BqT6YskXQC8rNfmI21f3ol4IkaDiUdd2ukQWnLvcW/sdAixHI3Jzsn221s9RtLPKI8FWZnyiJCP2P7bSMcWERGDy7TeUu+0vTWwBTAB2LfD8URErLDa1jlJWkvSpZLm1dV2+0naVtK1km6SdHlPyiBJH5B0Qy37E0lr1u371mPnSZpRt60u6XuSFki6RdJudfs0SedL+llNi/TFgeKz/af6cmVgVfq4CTfpiyIi2qOdI6c3APfb3tr2FsDPgG8A+9jeFjgN+Fwte77t7epI5nbgfXX7McBedftb6raPANjeEtgfOENSz2PipwD7AVsC+0lqvNH2OSRdTslG8QRwXu/9SV8UEdEe7bzmtAD4kqQvAJcAf6RMoV1Zc7WuBPTc+LqFpM9SspCPA3oWMswCTpf0Y+D8um1nSieH7Tsk/RbYrO67yvbjAJJ+BWzMsjn0lmF7r9qxnQ3sDlw53EZHjBZZYBDdpG2dk+07JW0L7A18nvKL/zbbO/ZR/HTgbbbnSZoG7FrrOETSDsAbgbmSptB3XrweLac2sv2UpIsomcrTOUVEdEA7rzltAPzF9lnAl4AdgAmSdqz7V5G0eS2+NvBAzUR+YEMdk2zPsX0M8AglH96MnjKSNgM2An7dYmzjGq53rUzpQO8YcmMjImJY2jmttyVwvKQlwNPAh4FngK9LWqfG8jVKlvL/AOZQHgq4gKX59I6XtClltHQVMI/SiZwkaUGtb5rtxQM81qkvawEXSVqNMr14NXDSMNoaERHDkPRFQ5T0RRERrUv6ooiIGLW6LkOEpGuAw233OyypiySm2v7oEOqfA6zWa/N7gFcA0yn3N82zfUCrdUeMZsszfVFWAkaruq5zWt5s79B7W72O9e/ATrb/KOnF7Y8sIiJ6DHtaT9InJB1WX39V0tX19esknSVpT0mzJd0s6VxJ4+r+PrNDNNT7PEln1PudkHSwpDslXQvs1FDuzZLm1OwQP5e0Xj32LkkTGur6jaR1+2nGB4Bv2f4jgO2Hhvu5RETE0I3ENacZwC719VRgXF0CvjNlpd2ngD1svwq4Efh43d9fdggoI7qzgTttf6p2XP9J6ZReD7yyoewvgFfb3gb4EfAJ20uAs1i6DH0PylTdI/20YTNgM0mzJF0n6Q19FUr6ooiI9hiJab2bgG0lrU256fVmSie1C+Xx56+kPL0WSs662cDf0392CIDvAD+23dNh7QBcY/thAEnnsDQLxIbAObUDWxVYWLefBvw3ZXn6PwHfG6ANKwObUm723RCYKWkL2481FrJ9MnAylNV6g380ERExFMPunGw/Lele4GDgl8B8YDdgEqWjuNL2/o3HSNqS/rNDUOvZTdKXbT/Vc6p+yn4D+IrtiyTtSlnUgO37JD0oaXdK53ZgP8cDLAKus/00sFDSrymd1Q0DHBMxpmTRQnSTkVpKPgM4vP47EzgEmAtcB+wkaRMASWvWLA6/pv/sEADfBS4Dzq0ZG+YAu0p6UZ0SbHycxTrA7+rrg3rFdSpleu/Hgzyb6UJKh0q9LrUZcE8L7Y+IiBE0Up3TTMqD+mbbfhB4CphZp+GmAT+UNJ/SWU22/VdgH+ALkuZROrL/11ih7a9QpgjPBB6kjIhmAz+v23tMp3RiMykpjRpdREkcO9CUHpTEso/W5LD/Axxh+9FmGx8RESNrTGeIkDQV+KrtXQYt3KJkiIiIaF2zGSLG7H1Oko6i5O8b6FpTRER0oXZmJb9M0vgWyk+UdOtQz2f7ONsb2/5FQ52flDS319cn+zj3n4d63oiIGL52Ps9p73ada4AYPsey91NFdK2/+5+5bT3f73eb0tbzRQxkxEZOTWSKuFfSunVEdLukUyTdJukKSWvUsttKmidpNvXx63X75pKuryOd+ZI2rfXcUbNIzJd0nqQ1G+p5TvYJSZMk/axunylpct3+sprF4gZJx47UZxIREUMzktN6A2WKmNmr7KaUdEGbA48B76jbvwcc1sf9T4cAJ9ieUuteVLf/PXCy7a2APwGHDpJ94mTgn+v2w4ET6/YTgG/b3g74/VA/gIiIGBkj2Tn1zhQxm6WZInp3Tgttz204bqLKAwfH2762bj+zofxs4GhJRwIb236ybr/P9qz6+ixKR9iYfWIuJX3ShjWn3/+jLDufS8lC0ZPPbyfgh32cdxlJXxQR0R4jds1pkEwRt/cqvrjh9d+ANShPt+1zXbvtH6g86uKNwOWS3k+5SbZ3edd6npN9QtLzgcfq6KvP0wzYQJK+KCKiXUZ6QURPpoh/oiR9/Qpwk21rkMem235M0uOSdq4r7J5dAi7p5cA9tr9eX29F6Zw2krSj7dnA/pQksM9mn7A9u07zbWb7NkkLJe1r+1yVgLayPQ+YBbyLZZPFRnRUFijEimykl5L3mSmiheMPBr5VF0Q82bB9P+DWOh03Gfh+3X47cFDNPvFCynWjgbJPHAi8r26/DXhr3f4vwEck3UBJhxQRER00ajNESJoIXGJ7i06cPxkiIiJa12yGiLbdhBsREdGsUZu+yPa9lFV5ERExxmTkFBERXWdUjZxqbr4DbJ9Y3+8KHG77TU0e/1HgY5Tl7RN6Htsu6UDgyFrsz8CH6yq+iI646upJbT/n63a/u+3njOjPaBs5jQcOHcbxs4A9gN/22r4QeG3NNHEs9V6miIjojLZ3Tg058U6VdKuksyXtIWmWpLskbS9puqTTJF0j6Z6enH3AccCkmmPv+LptXM2rd0etq98bqmzfUq9V9d7+S9t/rG+vAzbsJ/ZkiIiIaINOjZw2oeSz24py39IBlNRDhwNH1zKTgb2A7YFP15tpjwLutj3F9hG13DaUqbpXAi+npCIajvcBP+1rh+2TbU+1PXXChAnDPE1ERPSnU53TQtsLbC+h3Ax7lcsNVwuAibXMpbaE53OXAAASSUlEQVQX1+tCDwHr9VPX9bYX1brmNhzfMkm7UTqnIwcrGxERy0+nFkQ05tZb0vB+CUtj6p1/r79Ymy03IElbAacC/2D70aHUETFSsjghVnSjbUHEE8DaI12ppI2A84H32L5zpOuPiIjWjKrOqY5oZtWFFMcPekAvkg6TtIiy4GG+pFPrrmOAFwEn1sUWyUsUEdFBoza3Xqclt15EROuSWy8iIkatMdM51funbq2vL6jTc41fe3U6xoiIaM6oSl/ULNtv73QMEc2aPn16p0MAuieOCBhDI6dqJUmnSLpN0hWS1qhZJqYCSFq3PkoeSdMkXSjp4vqE3I9K+rikWyRdJ+mFHW1JRMQKbKx1TpsC37K9OfAY8I5Bym9ByU6xPfA54C+2twFmA+/tXTjpiyIi2mOsdU4Lbc+tr29i8GwR/2P7CdsPA48DF9ftjZkqnpX0RRER7THWOqe+skU8w9J2rj5A+f4yVURERJutCL+A7wW2Ba4H9ulsKBHPlYUIEc811kZOffkS8GFJvwTW7XQwERExuGSIGKJkiIiIaF0yRERExKiVzikiIrrOirAgIqIrLDpqZqdDGNCGx+3S6RAinjVmRk6Sxks6tOH9rpIuaeH4mQ15+O6XdOHyiTQiIgYzZjonYDxw6KCl+mF7F9tTbE+hZIg4f8Qii4iIlnRV51Qzi98h6dT6QMGzJe0haZakuyRtL2m6pNNqzrx7JB1WDz8OmFRHPj0PIhwn6bxa59mS1EQMawO7A88ZOSV9UUREe3RV51RtApwAbAVMpuS+2xk4HDi6lpkM7EXJifdpSasARwF319HPEbXcNsDHgFcCLwd2auL8bweusv2n3juSvigioj26sXNaaHuB7SXAbZSOwiyb7+5S24ttPwI8BKzXT13X215U65rL4Ln2APYHfjicBkRExPB042q9ZvLd9ZVDb7C6BioHgKQXUUZjeR5UjLishotoXjeOnIbqCWDtYdaxL3CJ7adGIJ6IiBiiMdM52X4UmFUXUhw/6AF9exeZ0ouI6Ljk1hui5NaLiGhdcutFRMSo1XULIiRNpFz32aLJ8qfX8udJOhX4iu1f9SozDZhq+6OSLgBe1quaI21fPtzYI3r78n5v6nQITfu3c5pOqBKx3HVd5zQctt/fRJmsxIuI6HLdOq23kqRTJN0m6QpJa0iaIuk6SfMlXSDpBb0PqlkjptbXB0u6U9K1NNx8K+nNkuZIukXSzyWtJ+l5NQPFhFrmeZJ+IykPJ4yI6IBu7Zw2Bb5le3PgMeAdwPcp029bUW7I/XR/B0taH/hPSqf0ekqGiB6/AF5texvgR8An6k26ZwEH1jJ7APPqTb6N9SZ9UUREG3Rr57TQ9tz6+iZgEjDe9rV12xnAawY4fgfgGtsP2/4rcE7Dvg2ByyUtAI4ANq/bTwPeW1//E/C93pUmfVFERHt0a+fUO7PD+CHU0d8a+W8A37S9JfAhYHUA2/cBD0randK5/XQI54yIiBEwWhZEPA78UdIutmcC7wGuHaD8HOCEmo7oT5TMD/PqvnWA39XXB/U67lTK9N6Ztv82UsHHiisr4CKGZrR0TlA6kpMkrQncAxzcX0HbD0iaTnku0wPAzcBKdfd04FxJvwOuY9ll5RdRpvOeM6UXERHtkwwRDepKv6/aHjRDZzJERES0rtkMEaNp5LRcSToK+DBLV+xFRESHdOuCiLazfZztjW3/otOxRESs6DJyihjAtw65utMhtM1HTtq90yFEPGvMjJwkjZd0aMP7XSU1vVRK0nclzasZKM6TNG75RBoREYMZM50T5V6oQwct1b9/tb11zUDxv8BHRyasiIhoVVd1TpImSrpD0qn1oYFnS9pD0qya+257SdMlnVbz6N0j6bB6+HHAJElzGx42OK6Ogu6odam/c9v+U41BwBr0cRNv0hdFRLRHV3VO1SbACcBWwGTgAGBn4HDg6FpmMrAXsD3waUmrAEcBd9ueYvuIWm4b4GOU3HovpyEBbF8kfQ/4fa3/G733J31RRER7dGPntND2gpqM9TbgKpebsRYAE2uZS20vrolZHwLW66eu620vqnXNbTi+T7YPBjYAbgf2G3ZLIiJiSLpxtV5jXr0lDe+XsDTe3rn3+mtHs+WeZftvks6hJIVNpogVXFawRXRGN46chuoJYO2hHKhik57XwJuBO0YwtoiIaEE3jpyGxPajdeHErZSM4pe2cLiAMyQ9v76eR8kWERERHZDcekOU3HoREa1rNrfeWJrWi4iIMWLMTOtJGg8cYPvE+n5X4HDbb+pV7gKWfUwGwJGUp+BOBZ4Grgc+ZPvp5R13jKzbJ7+i0yGMWq+44/ZOhxDxrDHTObE0Q8SJAxWy/fa+tktaCXh3ffsD4P3At0cywIiIaE5XTet1OEPEZa4oI6cNl3NzIyKiH13VOVUdyxABUOt6D/CzPvYlfVFERBt0Y+fUsQwR1YnADNsze+9I+qKIiPboxmtOHcsQIenTwATgQ80GG90lF/UjxoZu7JyGasgZIgAkvZ8yVfi6OtKKiIgO6cZpvSGx/Sgwqy6kOH7QA57rJMr04Oy6qOKYkY0wIiKalQwRQ5QMERERrUuGiIiIGLXG0jWnpvSXIcL25Z2IJyIinmvMdE7Npi8aIEPERyn3RE0CJtRl6jFMW56xZadDiCYtOGhBp0OIeNZYmtbrSV80VLOAPYDfjkw4ERExVF3VOXU4fdEttu9d3m2MiIjBdeO03ibAvsAHgRtYmr7oLZT0RXMp6Yt2o9zX9GtJ36akL9rC9hR4dlpvG2Bz4H7KyGgn4BdDDUzSB2tcbLTRRkOtJiIiBtFVI6eq0+mL+pX0RRER7dGNI6eOpS+KkZeL7BExFN04chqqYaUvioiI7jFmOqfhpi+SdJikRZTnOM2XdOqIBxkREU1J+qIhSvqiiIjWJX1RRESMWivcAoGkL4qI6H4rXOc0QPqizwHvBV5ge1x7oxrlpq/T6QhiJEx/vNMRRDwr03pLXQxs3+kgIiKijZ2TpLUkXSppXl1Rt5+kbSVdK+kmSZdLWr+W/YCkG2rZn0has27ftx47T9KMum11Sd+TtEDSLZJ2q9unSTpf0s9q6qMvDhSf7etsPzBIGz4o6UZJNz788MMj88FERMRztHPk9Abgfttb294C+BnwDWAf29sCpwGfq2XPt72d7a2B24H31e3HAHvV7W+p2z4CYHtLYH/gDEmr131TgP2ALYH9JL10OA1IhoiIiPZo5zWnBcCXJH0BuAT4I7AFcGXNx7oS0DNy2ULSZymZxscBPYsVZgGnS/oxcH7dtjOlk8P2HZJ+C2xW911l+3EASb8CNgbuW24tjIiIEdG2zsn2nZK2BfYGPg9cCdxme8c+ip8OvM32PEnTgF1rHYdI2gF4IzBX0hSg30zjJH1Re+RCekSMsHZec9oA+Ivts4AvATsAEyTtWPevImnzWnxt4AFJqwAHNtQxyfYc28cAjwAvBWb0lJG0GbAR8Os2NSsiIpaDdo4ktgSOl7QEeBr4MPAM8HVJ69RYvkbJRP4fwBzKg/8WsDRn3vGSNqWMlq4C5gF3ACdJWlDrm2Z78QCPbupTXTBxALBmTWN0qu3pQ29uREQMVdIXDVHSF0VEtC7piyIiYtRa4RYISJoDrNZr83ts58FDERFdou2dU119d4Xt++v7e4Gp9am2I3meyyjXkAAOsH0igO0d+ii7saSbKMvZVwG+YfukkYynt4lHXbo8q49o2b3HvbHTIUQ8qxPTetOADUaiIkn9dq6297b9GOVeqUMHqeoB4P/ZnkJZRXhUXV0YEREdMGjnJOkTkg6rr78q6er6+nWSzpK0p6TZkm6WdK6kcXX/MTUF0a2STlaxDzAVOFvSXElr1NP8cz1+gaTJ9fi1JJ1W67hF0lvr9mn1PBcDV0haX9KMWt+tknap5e6VtC5wHDCp7u/zIYS2/2q7556o1fr7XJK+KCKiPZoZOc0AdqmvpwLj6v1HO1OWeX8K2MP2q4AbgY/Xst+sKYi2ANYA3mT7vFrmQNtTbD9Zyz5Sj/82cHjd9kngatvbAbtRlpGvVfftCBxke3fK1N3lddSzNTC3V/xHAXfX8x3RXyMlvVTSfEoGiS/0TDs2SvqiiIj2aKZzugnYVtLalIwLsymd1C7Ak8ArKY9HnwscREkRBLCbpDn1/qPdgc2fU/NSPamIbgIm1td7UqbX5gLXAKtTbrAFuNL2H+rrG4CDJU0HtrT9RBNteg7b99neCtgEOEjSekOpJyIihm/QBRG2n66LFg4GfgnMp4xkJgELKR3F/o3H1MSrJ1IWOtxXO47V6V/PlFpjiiEB77C9TLaHmr7o/xrimyHpNZSURmdKOt729wdrV39s3y/pNkrne95Q6xlMLj5HRPSv2QURMyjTbTOAmcAhlOmz64CdJG0CIGnNmkKopyN6pF6D2qehridYmvFhIJdTrkWp1r1NX4UkbQw8ZPsU4LvAq3oVGfR8kjbsuf4l6QXATiQFUkRExzTbOc0E1gdm234QeAqYafthyuq7H9brNdcBk+squVMo16QupEy99Tidkm6ocUFEX46lLOueL+nW+r4vu1KSwN4CvAM4oXGn7Ucp04639rcgAngFMEfSPOBa4Eu57ykionOSvmiIJD1Myf3XaetSkuCOZmlDd0gbusNYb8PGtgddUZbOaZSTdGMzeaq6WdrQHdKG7pA2FCtU+iJJWwJn9tq8uK+sERER0TkrVOdUryNN6XQcERExsGQlH/1O7nQAIyBt6A5pQ3dIG8g1p4iI6EIZOUVERNdJ5xQREV0nndMoI+mFkq6UdFf99wX9lNtI0hWSbpf0K0kT2xtp/5ptQy37fEm/k/TNdsY4mGbaIGlKzdh/m6T5kvbrRKy9SXqDpF9L+o2ko/rYv5qkc+r+Od30s9OjiTZ8vP7cz5d0Vc0k01UGa0NDuX0kWVLXLS9vpg2S3lm/F7dJ+kHTldvO1yj6Ar4IHFVfH0XJoN5XuWuA19fX44A1Ox17q22o+08AfkDJct/x2FtpA7AZsGl9vQHluWHjOxz3SsDdwMuBVYF5wCt7lTkUOKm+fhdwTqc/7yG0Ybeen3ngw6OxDbXc2pS0cddRcpV2PPYWvw+bArcAL6jvX9xs/Rk5jT5vBc6or88A3ta7gKRXAivbvhLA9p9t/6V9IQ5q0DYASNoWWA+4ok1xtWLQNti+0/Zd9fX9wENAp5+1sj3wG9v32P4r8CNKWxo1tu084HU9OS67xKBtsP0/DT/z1wEbtjnGwTTzfYCStu2LlJRx3aaZNnwA+JbtPwLYfqjZytM5jT7r2X4AoP774j7KbAY8Jul8lQc1Hi9ppbZGObBB2yDpecCXgX6fwdVhzXwfniVpe8pfl3e3IbaBvITyzLIei+q2PsvYfgZ4HHhRW6JrTjNtaPQ+4KfLNaLWDdqGmuz6pbYvaWdgLWjm+7AZsJmkWZKuk/SGZitfoW7CHS0k/Rz4uz52fbLJKlamPPJjG+B/gXMoCXq/OxLxNWME2nAocJnLI1dGLrAWjEAbeupZn5KZ5CDbS0YitmHo68PsfT9JM2U6qen4JL2b8vy51y7XiFo3YBvqH2dfpfy/7VbNfB9Wpkzt7UoZvc6UtIVLcvABpXPqQrb36G+fpAclrW/7gfpLr69h8iLgFtv31GMuBF5NGzunEWjDjsAukg6lXDNbVdKfbfd74XikjUAbkPR84FLgU7avW06htmIR8NKG9xsCvZ/63FNmkaSVgXWAP9A9mmkDkvag/CHxWtuLe+/vsMHasDawBXBN/ePs74CLJL3F9o1ti3Jgzf4sXWf7aWChpF9TOqsbGESm9UafiyhPHKb++999lLkBeIGknusbuwO/akNszRq0DbYPtL2R7YmUZ4l9v50dUxMGbYOkVYELKLGf28bYBnIDsKmkl9X43kVpS6PGtu0DXO16NbtLDNqGOiX2HeAtrVznaKMB22D7cdvr2p5Y/w9cR2lLt3RM0NzP0oWUxSlIWpcyzXdPU7V3esVHvlpeIfMi4CrgrvrvC+v2qcCpDeVeT3lq8QLKM7RW7XTsrbahofw0um+13qBtAN4NPE15MGfP15QuiH1v4E7K9a9P1m2fofzyg/Kw0HOB3wDXAy/vdMxDaMPPgQcbPveLOh1zq23oVfYaumy1XpPfBwFfofxxvAB4V7N1J31RRER0nUzrRURE10nnFBERXSedU0REdJ10ThER0XXSOUVERNdJ5xQREV0nnVNERHSd/w9+QmRxGyQ7RwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21213787c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl0nfV95/H3R5LlfcMbeJFswCwmGGOEwUBCU0wCJDFZmrDbaelw0hzOtMPJTOkwh5wh0zMNND0zndIZaCcTG0jYAqlLTAhxCEljGbxisM1ijLXYxgveLdvavvPHfWQuQvaVsB7de6XP6xwd3fs8v+c+Xz0YffSsX0UEZmZmJ1KS7wLMzKzwOSzMzCwnh4WZmeXksDAzs5wcFmZmlpPDwszMcnJYmJlZTg4LMzPLyWFhZmY5leW7gO4yevTomDx5cr7LMDMrKitXrtwVEWNyjes1YTF58mRWrFiR7zLMzIqKpJrOjPNhKDMzy8lhYWZmOTkszMwsJ4eFmZnl5LAwM7OcHBZmZpaTw8LMzHJyWJiZFbFnV9fz1Io60m6RnWpYSLpG0luSNkq6u4P5d0laL2mtpCWSKrPmVUj6paQNyZjJadZqZlZsmlta+dsX3ubZ1VuQlOq6UgsLSaXAg8C1wDTgJknT2g1bDVRFxHTgaeD+rHkLgQci4lxgFrAjrVrNzIrRkjd3sGXvYebNnpz6utLcs5gFbIyITRHRCDwOXJ89ICJeioiG5O0yYCJAEiplEfFiMu5g1jgzMwMWVm9m/PABzDl3bOrrSjMsJgB1We/rk2nHczvwfPL6LGCvpGckrZb0QLKnYmZmwMYdB/j9xg+45dJKykrTP/2c5ho6OoDW4RkYSbcCVcADyaQy4NPAd4CLgdOBb3aw3B2SVkhasXPnzu6o2cysKCysrqG8tIQbL57UI+tLMyzqgeyfYiKwtf0gSXOAe4C5EXE0a9nVySGsZuBnwMz2y0bEwxFRFRFVY8bkfMKumVmvcOBIEz9dWc8XLziNUUP698g60wyL5cBUSVMklQM3AouyB0i6EHiITFDsaLfsSEltCfCHwPoUazUzKxrPrt7CocYW5vfAie02qYVFskdwJ/ACsAF4MiLWSbpP0txk2APAEOApSWskLUqWbSFzCGqJpNfJHNL6p7RqNTMrFhHBgqWbuWDicC6YNKLH1ptq86OIWAwsbjft3qzXc06w7IvA9PSqMzMrPkvf/YB3dx7iB1+/oEfX6zu4zcyKyIKlmzllcDlfmH5aj67XYWFmViTq9zTwqw3bufHiSQzo17N3EzgszMyKxGOv1AJwy6WVOUZ2P4eFmVkRONLUwuOv1nL1tHFMGDGwx9fvsDAzKwLPrd3GnoamHr1cNpvDwsysCCys3syZY4cw+4xReVm/w8LMrMCtqdvL2vp9zJ9dmfqjyI/HYWFmVuAWLt3MkP5lfGXmxLzV4LAwMytguw4e5bm12/jazAkM6Z/qfdQn5LAwMytgTyyvo7GlldvydGK7jcPCzKxANbe08uiyGq44czRnjh2S11ocFmZmBepXG7azbd8R5s3u+Zvw2nNYmJkVqAVLa5gwYiBXnTsu36U4LMzMCtHb2w9QvekDbr20ktKS/Fwum81hYWZWgBZWb6a8rIQbeqhtai4OCzOzArP/SBPPrNrCl6aP55TB5fkuB3BYmJkVnGdW1tPQ2ML8y/J/YruNw8LMrIC0tgYLq2uYMWkE0yf2XNvUXBwWZmYF5Pfv7mLTrkMFtVcBKYeFpGskvSVpo6S7O5h/l6T1ktZKWiKpMmtei6Q1ydeiNOs0MysUC5bWMGpwOded37NtU3NJLSwklQIPAtcC04CbJE1rN2w1UBUR04Gngfuz5h2OiBnJ19y06jQzKxR1uxtY8uZ2bppVQf+ynm2bmkuaexazgI0RsSkiGoHHgeuzB0TESxHRkLxdBuTvkYpmZnn26Cs1lEjcfElFvkv5mDTDYgJQl/W+Ppl2PLcDz2e9HyBphaRlkr7c0QKS7kjGrNi5c+fJV2xmlidHmlp4Ynkdn5s2jvF5aJuaS5rPu+3olsPocKB0K1AFXJk1uSIitko6Hfi1pNcj4t2PfFjEw8DDAFVVVR1+tplZMVj02lb2NjQxL89Plz2eNPcs6oHsWw8nAlvbD5I0B7gHmBsRR9umR8TW5Psm4DfAhSnWamaWNxHBwurNnDVuCJeefkq+y+lQmmGxHJgqaYqkcuBG4CNXNUm6EHiITFDsyJo+UlL/5PVo4HJgfYq1mpnlzeq6vbyxZT+3zZ6ct7apuaR2GCoimiXdCbwAlAI/jIh1ku4DVkTEIuABYAjwVLKBapMrn84FHpLUSibQ/iYiHBZm1istXLqZof3L+OqFJzqtm1+p9uiLiMXA4nbT7s16Pec4yy0Fzk+zNjOzQrDzwFF+/vo2brmkksF5bJuai+/gNjPLo8dfraWpJbitABocnYjDwswsT5paWnnslVo+PXU0Z4zJb9vUXBwWZmZ58uL67by//wjzC/Ry2WwOCzOzPFmwdDMTRw7ks+eMzXcpOTkszMzy4M339/PKe7sLpm1qLg4LM7M8eKS6hv5lJdxQVRhtU3NxWJiZ9bB9hzNtU+deMJ6RBdI2NReHhZlZD/vpynoON7Uw/7LJ+S6l0xwWZmY9qLU1eGRZDTMrRvCpCcPzXU6nOSzMzHrQ7zbu4r1dh4pqrwIcFmZmPWrh0s2MHtKfaz9VWG1Tc3FYmJn1kNoPGvj1Wzu4edYkysuK69dvcVVrZlbEPmybWtjPgeqIw8LMrAccbsy0Tf38eeM4dfiAfJfTZQ4LM7Me8K+vbWXf4cJtm5qLw8LMLGURwY+WbubscUO5ZEphtk3NxWFhZpayVbV7WL9tP/MuqyzYtqm5OCzMzFK2YGkNQweU8eUZhds2NZdUw0LSNZLekrRR0t0dzL9L0npJayUtkVTZbv4wSVsk/UOadZqZpWXH/iMsfn0bX79oUkG3Tc0ltbCQVAo8CFwLTANukjSt3bDVQFVETAeeBu5vN/97wMtp1WhmlrafvFpHc2vht03NJc09i1nAxojYFBGNwOPA9dkDIuKliGhI3i4DJrbNk3QRMA74ZYo1mpmlJtM2tYYrzxrDlNGD813OSUkzLCYAdVnv65Npx3M78DyApBLgB8B/TK06M7OUvbDufXYcOMr8y4p7rwIgzQNoHZ3yjw4HSrcCVcCVyaRvA4sjou5EVw5IugO4A6CiouKkijUz624Ll9Yw6ZSBXHlW4bdNzSXNsKgHsltATQS2th8kaQ5wD3BlRBxNJs8GPi3p28AQoFzSwYj4yEnyiHgYeBigqqqqwyAyM8uHDdv28+rm3fzn684pirapuaQZFsuBqZKmAFuAG4GbswdIuhB4CLgmIna0TY+IW7LGfJPMSfCPXU1lZlaoFiZtU79RJG1Tc0ntnEVENAN3Ai8AG4AnI2KdpPskzU2GPUBmz+EpSWskLUqrHjOznrKvoYmfrd7Cl2dMYMSg4mibmkuqF/1GxGJgcbtp92a9ntOJz/gR8KPurs3MLC1PrazjcFNL0V8um813cJuZdaO2tqlVlSOLqm1qLg4LM7Nu9PI7O6n5oIF5RdY2NReHhZlZN1q4dDNjhvbnmvNOzXcp3cphYWbWTTbvOsRv3t7JTbMqiq5tai6966cxM8ujR5fVUCpxyyW97yZhh4WZWTc43NjCkyvq+PynTmXcsOJrm5qLw8LMrBv8y5ot7D/SzPwibZuai8PCzOwkRQQLqms459ShXDx5ZL7LSYXDwszsJK2o2cOGbfuZf9nkom2bmovDwszsJC1YuplhA8q4fsb4fJeSGoeFmdlJ2L7/CL94432+UTWJQeXF2zY1F4eFmdlJ+PErtbREcOulvec5UB1xWJiZfUKNza38+NVarjxrDJOLvG1qLg4LM7NP6Bfr3mfngaO99nLZbA4LM7NP6JHqzVScMogrzxqT71JS57AwM/sE1m3dx/LNe5g3u5KSXtA2NReHhZnZJ/BIdQ0D+pXw9Yt6R9vUXBwWZmZdtLehkZ+t2cJXLpzA8EH98l1Oj3BYmJl10VMr6jnS1Mptl07Odyk9JtWwkHSNpLckbZR0dwfz75K0XtJaSUskVSbTKyWtlLRG0jpJ30qzTjOzzmpJ2qbOmnwK08YPy3c5PSa1sJBUCjwIXAtMA26SNK3dsNVAVURMB54G7k+mbwMui4gZwCXA3ZJ67330ZlY0Xn57B7W7G5h3We++Ca+9NPcsZgEbI2JTRDQCjwPXZw+IiJcioiF5uwyYmExvjIijyfT+KddpZtZpC5bWMHZofz7fy9qm5pLmL+EJQF3W+/pk2vHcDjzf9kbSJElrk8/4fkRsbb+ApDskrZC0YufOnd1UtplZx97bdYiX397JzZdU0K+0b/0Nm+ZP29GFx9HhQOlWoAp44NjAiLrk8NSZwHxJ4z72YREPR0RVRFSNGdP7b4oxs/x6dFkNZSXi5lm9r21qLmmGRT2QfQHyRKCjvYM5wD3A3KxDT8ckexTrgE+nVKeZWU4Njc08uaKOa88/jbG9sG1qLmmGxXJgqqQpksqBG4FF2QMkXQg8RCYodmRNnyhpYPJ6JHA58FaKtZqZndDPVm/lwJFm5s/uWye226T28PWIaJZ0J/ACUAr8MCLWSboPWBERi8gcdhoCPJV0l6qNiLnAucAPJAWZw1l/GxGvp1WrmdmJRAQLqzcz7bRhXFTZO9um5tLpsJB0BTA1Iv6fpDHAkIh470TLRMRiYHG7afdmvZ5znOVeBKZ3tjYzszS9+t5u3nz/AN//2vm9tm1qLp06DCXpu8BfAn+VTOoHPJpWUWZmhWRhdQ3DB/Zj7gUnuqCzd+vsOYuvAHOBQ3DspPPQtIoyMysU7+87wi/Wvc83qiYysLw03+XkTWfDojEiguTSV0m9uyWUmVnix6/U0NoH2qbm0tmweFLSQ8AISf8O+BXwT+mVZWaWf5m2qXV89uyxVI7q238jd+oEd0T8raSrgf3A2cC9yUloM7Ne6/k3trHr4FHm9dHLZbN1KiySw06/jogXJZ0NnC2pX0Q0pVuemVn+LKyuYfKoQXxmqp8Q0dnDUL8F+kuaQOYQ1B8DP0qrKDOzfHtjyz5W1uzhttmT+0Tb1Fw6GxZKng77VeB/RcRXyDx23MysV1pYvZmB/Ur5o4sm5ruUgtDpsJA0G7gF+HkyLbW7v83M8mnPoUb+Zc1WvjJzAsMH9o22qbl0Niz+HLgbeCZ5ZMcU4NfplWVmlj9PrqjjaHOrT2xn6ezeQQPQSqbb3a1kntfU4ePGzcyK2bG2qVNO4ZxT+07b1Fw6GxaPAd8B3iATGmZmvdJLb+6gfs9h/urac/NdSkHpbFjsjIh/TbUSM7MCsKB6M+OG9edz532s31qf1tmw+K6kfwaWAMcaFEXEM6lUZWaWB5t2HuR37+zirqvP6nNtU3PpbFj8MXAOmafNth2GCsBhYWa9xiPLauhXKm6cNSn34D6ms2FxQUScn2olZmZ5dOhoM0+vqOe6809j7NC+1zY1l87uZy2T5JvwzKzXenb1Fg4cbWbe7Mn5LqUgdXbP4gpgvqT3yJyzEBAR4W52Zlb02tqmfmrCMGZWjMh3OQWps3sW1wBTgc8BXwK+mHw/IUnXSHpL0kZJd3cw/y5J6yWtlbREUmUyfYakaknrknk3dP5HMjPrmmWbdvP29oPMu3Ryn22bmktnH1Fe09UPllQKPAhcDdQDyyUtioj1WcNWA1UR0SDpz4D7gRvI3AQ4LyLekTQeWCnphYjY29U6zMxyWVi9mRGD+jF3xvh8l1Kw0rw2bBawMSI2RUQj8DhwffaAiHgpeUAhwDJgYjL97Yh4J3m9FdgB+BnBZtbttu49zC/Xb+eGqkkM6Nd326bmkmZYTADqst7XJ9OO53bg+fYTJc0CyoF3u7U6MzPgx6/Uum1qJ6T55NiODvx1+Dyp5HlTVcCV7aafBjwCzI+Ijz1mRNIdwB0AFRUVJ1uvmfUxR5tbeHx5LVedM5ZJpwzKdzkFLc09i3og+86WicDW9oMkzQHuAeZGxNGs6cPIPA79v0TEso5WEBEPR0RVRFSNGeOjVGbWNc+//j67Djb6ctlOSDMslgNTJU2RVA7cCCzKHiDpQuAhMkGxI2t6OfAssDAinkqxRjPrwxZUb+b00YO54szR+S6l4KUWFhHRDNwJvABsAJ5MemHcJ2luMuwBYAjwlKQ1ktrC5BvAZ4BvJtPXSJqRVq1m1vesrd/L6tq93Da70m1TOyHVbncRsRhY3G7avVmv5xxnuUeBR9Oszcz6toXVNQwqL+VrbpvaKX6sopn1ObsPNbLota185cIJDBvgtqmd4bAwsz7nieV1NDa3+sR2FzgszKxPaWkNHl1Ww6Wnn8LZpw7NdzlFw2FhZn3Kkg3b2bL3MPO9V9ElDgsz61MeWVbDacMHcPU0t03tCoeFmfUZG3dk2qbeckkFZW6b2iXeWmbWZzy6rIby0hJunOXHA3WVw8LM+oSDR5t5emU9X5h+GqOH9M93OUXHYWFmfcKzq+o5eLSZebP9dNlPwmFhZr1eRLCguobzJwxnxiS3Tf0kHBZm1utVv/sBG3ccZN7sSrdN/YQcFmbW6y2o3szIQf340gVum/pJOSzMrFfbsvcwL67fzg0XV7ht6klwWJhZr/bjV2oAuOUSXy57MhwWZtZrHWlq4Sev1nHVuePcNvUkOSzMrNda/Po2dh9q9HOguoHDwsx6rQXVNZw+ZjCXnzkq36UUPYeFmfVKa+r28lrdXuZd6stlu4PDwsx6pYXVmxnstqndJtWwkHSNpLckbZR0dwfz75K0XtJaSUskVWbN+4WkvZKeS7NGM+t9Pjh4lOde28ZXZ05kqNumdovUwkJSKfAgcC0wDbhJ0rR2w1YDVRExHXgauD9r3gPAbWnVZ2a91+PL62hsafVzoLpRmnsWs4CNEbEpIhqBx4HrswdExEsR0ZC8XQZMzJq3BDiQYn1m1gs1t7Ty2LIaLjtjFFPHuW1qd0kzLCYAdVnv65Npx3M78HyK9ZhZH7DkzR1s3XeEeb5ctluVpfjZHV1+EB0OlG4FqoAru7QC6Q7gDoCKCt+daWaZE9vjhw9gzrlj811Kr5LmnkU9MCnr/URga/tBkuYA9wBzI+JoV1YQEQ9HRFVEVI0ZM+akijWz4rdxxwF+v/EDbrm00m1Tu1maW3M5MFXSFEnlwI3AouwBki4EHiITFDtSrMXM+oCF1Unb1Isn5R5sXZJaWEREM3An8AKwAXgyItZJuk/S3GTYA8AQ4ClJayQdCxNJvwOeAq6SVC/p82nVambF78CRJn66sp4vTj+NUW6b2u3SPGdBRCwGFrebdm/W6zknWPbTKZZmZr3MM6u2cKixhXmXTc53Kb2SD+qZWdHLtE3dzAUT3TY1LQ4LMyt6v9/4AZt2HvLlsilyWJhZ0VtQvZlTBpfzhemn5buUXsthYWZFrX5PA0s2bOfGiye5bWqKHBZmVtQee6UWgFsu9XOg0uSwMLOidaSphcdfreXqaeOYMGJgvsvp1RwWZla0nlu7jT0NTW6b2gMcFmZWlCKCBUs3c+bYIcw+w21T0+awMLOitKZuL69v2ce82W6b2hMcFmZWlBZW1zCkfxlfnem2qT3BYWFmRWfngaP8fO02vjZzAkP6p/rUIks4LMys6DyxvJbGllZu84ntHuOwMLOi0tzSymOv1HLFmaM5c+yQfJfTZzgszKyo/GrDdrbtO8K82b4Jryc5LMysqCxYWsOEEQO56txx+S6lT3FYmFnReHv7Aao3fcAtl1ZQWuLLZXuSw8LMisbC6s2Ul5VwQ5XbpvY0h4WZFYX9R5p4ZtUWvjR9vNum5oHDwsyKwk9X1tPQ2ML8y3xiOx9SDQtJ10h6S9JGSXd3MP8uSeslrZW0RFJl1rz5kt5JvuanWaeZFbbW1uCR6hpmTBrB9Ilum5oPqYWFpFLgQeBaYBpwk6Rp7YatBqoiYjrwNHB/suwpwHeBS4BZwHcljUyrVjMrbP+2cRebdh3yXkUepblnMQvYGBGbIqIReBy4PntARLwUEQ3J22VA20NePg+8GBG7I2IP8CJwTYq1mlkBW1hdw6jB5Vx3vtum5kuaYTEBqMt6X59MO57bgee7sqykOyStkLRi586dJ1mumRWiut0NLHlzOzfNqqB/mdum5kuaYdHRRdDR4UDpVqAKeKAry0bEwxFRFRFVY8aM+cSFmlnhevSVGkokbr6kIt+l9GlphkU9kH0x9ERga/tBkuYA9wBzI+JoV5Y1s97tSFMLTyyv4+pzxzHebVPzKs2wWA5MlTRFUjlwI7Aoe4CkC4GHyATFjqxZLwCfkzQyObH9uWSamfUhi17byt6GJub5xHbepfYg+IholnQnmV/ypcAPI2KdpPuAFRGxiMxhpyHAU0mnq9qImBsRuyV9j0zgANwXEbvTqtXMCk9b29SpY4cw+3S3Tc23VLuGRMRiYHG7afdmvZ5zgmV/CPwwverMrJCtqt3Luq37+d6XP+W2qQXAd3CbWUFaWL2Zof3L+OqFJ7qI0nqKw8LMCs7OA0dZ/Po2vnbRRAa7bWpBcFiYWcF5/NVamlqC29zgqGA4LMysoDQlbVM/PXU0Z4xx29RC4bAws4Ly4vrtvL//CPNnT853KZbFBwPNLK+aWlp5c9sBVtXuYWXNHv5t4y4mjBjIZ88Zm+/SLIvDwsx61O5Djayq2XMsHNbW7+NwUwsA44b159LTT+H2K6a4bWqBcViYWWpaWoN3dhxgVc1eVtbsYXXtHjbtOgRAWYk4b/wwbrh4EhdVjmRm5UjGDx/geyoKlMPCzLrN/iNNrKnNBMOq2j2sqd3LgaPNAIwaXM6FFSP5etUkZlZkmhgNLPdTZIuFw8LMPpGI4L1dh5Jg2Muqmj28veMAEVAiOGvcUObOGM/MipFcVDmSylGDvNdQxBwWZtYpDY3NvFa3j1W1e46dc9jT0ATA0AFlzKwYyRemn8bMipFcMGk4Qwf0y3PF1p0cFmb2MRFB/Z7Dx4JhZe0eNmw7QEtrpq3MGWMGc/W0ccf2Gs4YM4QSn5Du1RwWZsaRphbWbd137ET0qto97DiQaS8zqLyUGZNG8GdXnsFFlSO5sGIEIwaV57li62kOC7M+aPv+I5k9hiQY3tiyn8aWVgAqThnEZWeMSoJhJOecOpSyUt+/29c5LMx6ubab3lbW7GZVcqXSlr2HASgvK2H6hOH88eWTmZnsNYwdOiDPFVsh6vNhERG0tAalJfKVGtYrtN30tjI53/Ba/V6ONGX2Gk4dNoCLKkfyJ1dMYWbFCM4bP5zyMu81WG59Pix2H2rkov/2KyBzuV9ZSQmlJfr4lzLfy0o/fH2iMSUSZR8bU0JZyYfzSko+OiZ7WqfHHFtXyQnHZNeXa0xZSQklJXz0u3CYFqC2m95W1uxhVc1eVtXu4b12N73dNKvi2Ilo97G2T6rPh8WAfqXcdfVZtLRm9jCaW4PWCJpbku+trbS0Qktra2Zeh2PiI8sfbWo94ZiPr6s1My0+nJdcdFJQsoNmUHkZIwb1Y/jAHF+D+jEieT1sYD8G9PNNWCdj3+Em1tTtPXbpavub3mZWjuQbVZk7os+fMNw3vVm3STUsJF0D/E8yPbj/OSL+pt38zwD/A5gO3BgRT2fN+z7wheTt9yLiiTRqHNy/jH9/1dQ0PvqktB0eOxYyEbS0fDRQskPnI9MiaElCLhN2H1+m7XOaj/OZHY9JPrOllUONzew73MTehia27z/C29sPsO9wEweONJ/w5+pfVvKxkBk2sB8jBpYn78sYfmx++UfG9bXDJRHBpl2HjgXDqpq9H7np7exThzF3xnguqszsNVSc4pveLD2phYWkUuBB4GqgHlguaVFErM8aVgt8E/hOu2W/AMwEZgD9gZclPR8R+9Oqt9AoOTRUVmR/GLa0BvsPN7Gv3dfew00fTm9om9bIlr1H2LAtEzQHj544aAb2Kz0WNMOSABnRbi/meHs4xXA1T/ZNb23PUWq76W3YgDIuTG56u6hyJBdMGsEQd5CzHpTmv7ZZwMaI2AQg6XHgeuBYWETE5mRea7tlpwEvR0Qz0CzpNeAa4MkU67VuUFoiRg4uZ+Tgrl+H39TS2mHQZAdMW/DsO9xE3e4G3kheNzS2nPCzh/QvywqZsiRoyo8FTIfhk0xP4+mn2Te9tV2+2tFNbxdVjmRmhW96s/xLMywmAHVZ7+uBSzq57GvAdyX9HTAI+CxZIWO9U7/SEkYN6c+oIf27vGxjc+tHwmV/sueSCZnmY3sybWH03q5D7Du8l70NTRxtbv+3ykcNHVD2kQBpv3dzLHjaBc3QAWXHfsG33fTWdiJ6Ze0edra76e3bf3AGMyt805sVpjTDoqM/gzp12jYifinpYmApsBOoBj52jELSHcAdABUVFZ+8Uit65WUljBnanzFDux40R5paPrJHs7fhOHs3ydfb2w8e29tpu5GtIxIMG9CPYQPL2L7v6EduervizNHMrBjBzMqRnD3ON71Z4UszLOqBSVnvJwJbO7twRPw18NcAkn4MvNPBmIeBhwGqqqoK8PohKwYD+pUyoF8pY4d17Wa0iOBIU2tWyDR2sHeT+T5u2ABmJoeUPkmgmeVbmmGxHJgqaQqwBbgRuLkzCyYnx0dExAeSppO5WuqXqVVq9glIYmB5KQPLSzl1uO96tt4ttbCIiGZJdwIvkLl09ocRsU7SfcCKiFiUHGp6FhgJfEnSf42I84B+wO+SywD3A7cmJ7vNzCwPUr32LiIWA4vbTbs36/VyMoen2i93hMwVUWZmVgB8Vs3MzHJyWJiZWU4OCzMzy8lhYWZmOTkszMwsJ4eFmZnlpIjeceOzpJ1AzUl8xGhgVzeV051cV9e4rq5xXV3TG+uqjIgxuQb1mrA4WZJWRERVvutoz3V1jevqGtfVNX25Lh+GMjOznBwWZmaWk8PiQw/nu4DjcF1d47q6xnV1TZ+ty+cszMwsJ+9ZmJlZTn02LCQ9IOlNSWslPStpxHHGXSPpLUkbJd3dA3V9XdI6Sa2Sjnt1g6TNkl6XtEbSigKqq6e31ymSXpT0TvJ95HHGtSTbao2kRSkuUEfoAAAGyklEQVTWc8KfX1J/SU8k81+RNDmtWrpY1zcl7czaRn/aAzX9UNIOSW8cZ74k/X1S81pJM9OuqZN1/YGkfVnb6t6OxqVQ1yRJL0nakPy/+OcdjElvm0VEn/wCPgeUJa+/D3y/gzGlwLvA6UA5md7g01Ku61zgbOA3QNUJxm0GRvfg9spZV5621/3A3cnruzv675jMO9gD2yjnzw98G/g/yesbgScKpK5vAv/QU/+eknV+BpgJvHGc+dcBz5Np0Xwp8EqB1PUHwHM9ua2S9Z4GzExeDwXe7uC/Y2rbrM/uWUTEL+PDhkrL6KCvBjAL2BgRmyKiEXgcuD7lujZExFtpruOT6GRdPb69ks9fkLxeAHw55fWdSGd+/ux6nwauUtLlK8919biI+C2w+wRDrgcWRsYyYISk0wqgrryIiG0RsSp5fQDYAExoNyy1bdZnw6KdPyGTxu1NAOqy3tfz8f84+RLALyWtlHRHvotJ5GN7jYuIbZD5nwkYe5xxAyStkLRMUlqB0pmf/9iY5I+VfcColOrpSl0AX0sOXTwtaVLKNXVGIf//N1vSa5Kel3ReT688OXx5IfBKu1mpbbNUO+Xlm6RfAad2MOueiPiXZMw9QDPwWEcf0cG0k758rDN1dcLlEbFV0ljgRUlvJn8R5bOuHt9eXfiYimR7nQ78WtLrEfHuydbWTmd+/lS2UQ6dWee/Aj+JiKOSvkVm7+cPU64rl3xsq85YReYRGQclXQf8DJjaUyuXNAT4KfAXEbG//ewOFumWbdarwyIi5pxovqT5wBeBqyI54NdOPZD9F9ZEYGvadXXyM7Ym33dIepbMoYaTCotuqKvHt5ek7ZJOi4htye72juN8Rtv22iTpN2T+KuvusOjMz982pl5SGTCc9A955KwrIj7IevtPZM7j5Vsq/55OVvYv6IhYLOkfJY2OiNSfGSWpH5mgeCwinulgSGrbrM8ehpJ0DfCXwNyIaDjOsOXAVElTJJWTOSGZ2pU0nSVpsKShba/JnKzv8MqNHpaP7bUImJ+8ng98bA9I0khJ/ZPXo4HLgfUp1NKZnz+73j8Cfn2cP1R6tK52x7Xnkjkenm+LgHnJFT6XAvvaDjnmk6RT284zSZpF5vfoBydeqlvWK+D/Ahsi4u+OMyy9bdbTZ/QL5QvYSObY3prkq+0KlfHA4qxx15G56uBdModj0q7rK2T+OjgKbAdeaF8XmataXku+1hVKXXnaXqOAJcA7yfdTkulVwD8nry8DXk+21+vA7SnW87GfH7iPzB8lAAOAp5J/f68Cp6e9jTpZ139P/i29BrwEnNMDNf0E2AY0Jf+2bge+BXwrmS/gwaTm1znB1YE9XNedWdtqGXBZD9V1BZlDSmuzfm9d11PbzHdwm5lZTn32MJSZmXWew8LMzHJyWJiZWU4OCzMzy8lhYWZmOTksrM+TdPAkl386uTP8RGN+oxM8rbezY9qNHyPpF50db3YyHBZmJyF5LlBpRGzq6XVHxE5gm6TLe3rd1vc4LMwSyV2vD0h6Q5leITck00uSRzqsk/ScpMWS/ihZ7Bay7hqX9L+TBxauk/Rfj7Oeg5J+IGmVpCWSxmTN/rqkVyW9LenTyfjJkn6XjF8l6bKs8T9LajBLlcPC7ENfBWYAFwBzgAeSx2B8FZgMnA/8KTA7a5nLgZVZ7++JiCpgOnClpOkdrGcwsCoiZgIvA9/NmlcWEbOAv8iavgO4Ohl/A/D3WeNXAJ/u+o9q1jW9+kGCZl10BZknr7YA2yW9DFycTH8qIlqB9yW9lLXMacDOrPffSB4ZX5bMm0bm8QzZWoEnktePAtkPhGt7vZJMQAH0A/5B0gygBTgra/wOMo9cMUuVw8LsQ8drQnSi5kSHyTzvCUlTgO8AF0fEHkk/apuXQ/Yzd44m31v48P/P/0DmeVwXkDkacCRr/ICkBrNU+TCU2Yd+C9wgqTQ5j/AZMg/7+zcyjYFKJI0j01azzQbgzOT1MOAQsC8Zd+1x1lNC5omzADcnn38iw4FtyZ7NbWTapLY5i8J44rD1ct6zMPvQs2TOR7xG5q/9/xQR70v6KXAVmV/Kb5PpTrYvWebnZMLjVxHxmqTVZJ5Iugn4/XHWcwg4T9LK5HNuyFHXPwI/lfR1Mk+EPZQ177NJDWap8lNnzTpB0pDIdEYbRWZv4/IkSAaS+QV+eXKuozOfdTAihnRTXb8Fro+IPd3xeWbH4z0Ls855TtIIoBz4XkS8DxARhyV9l0yf49qeLCg5VPZ3DgrrCd6zMDOznHyC28zMcnJYmJlZTg4LMzPLyWFhZmY5OSzMzCwnh4WZmeX0/wGhAkhEuHkHcgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x212196a3e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot important coefficients\n",
    "#   ridge.coef_ 就是训练出的系数w\n",
    "\n",
    "coefs = pd.Series(ridge.coef_, index = x_train.columns)\n",
    "print(\"Ridge picked \" + str(sum(coefs != 0)) + \" features and eliminated the other \" +  \\\n",
    "      str(sum(coefs == 0)) + \" features\")\n",
    "\n",
    "#正系数值最大的10个特征和负系数值最小（绝对值大）的10个特征\n",
    "imp_coefs = pd.concat([coefs.sort_values().head(10),\n",
    "                     coefs.sort_values().tail(10)])\n",
    "imp_coefs.plot(kind = \"barh\")\n",
    "plt.title(\"Coefficients in the Ridge Model\")\n",
    "plt.show()\n",
    "\n",
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1)) \n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 Linear Regression with Lasso regularization (L1 penalty) Lasso(L1)正则"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 0.0008015891035234101\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8XNWd9/HPT11WtVVc5CI3uWIMyJjYBkNohk1gSSghhBAC6w27JLtJdlMedkMSdp9nk2yyyW4KIQkhLC0JLcYQMN2Y4oaNe5GbLNvqvUuj8/wxY0WAZI1tXd0Z6ft+veZlzcyZO7+ra8137rnnnmvOOURERABi/C5AREQih0JBRES6KRRERKSbQkFERLopFEREpJtCQUREuikURESkm0JBRES6KRRERKRbnN8FnKzs7GyXn5/vdxkiIlFl48aNlc65nP7aRV0o5Ofns2HDBr/LEBGJKmZ2KJx26j4SEZFuCgUREemmUBARkW4KBRER6aZQEBGRbgoFERHpplAQEZFuUXeegojIUNfaEaC2uYOa5nbKG9ooq2ulrL6V+RMzOX96v+efnRaFgojIINh2pI7vPLOd5vYAZuAcBLocgS5HZ5ejtSNAS0eA5vYA7Z1dvS7jjgunKhRERKJdS3uALz26ifrWDuZPyMS54OOxMUZcrBEbE0NyfAzJ8bEkJcSSkRxPZnICGcnx5KYnMiY9iZy0RJLiYz2vVaEgIuKx//jzTvZXNvHI7QtZNC3b73JOSAeaRUQ89MbeCn739iE+v3hyxAcCKBRERDxT19zBP/9xC9NyU/nashl+lxMWdR+JiHjkkXXFlNa3suKziwfleMBA0J6CiIhHVu+pYOaYNOaNz/S7lLB5Fgpmdr+ZlZvZthO0udDMNpvZdjN73ataREQGW0t7gI2HarigwNshpAPNyz2FB4BlfT1pZpnAz4GrnHNzgOs8rEVEZFCtPVBFe6CLJVFwcLknz0LBObcaqD5Bk08DTzrnikPty72qRURksL2xt5KEuBjOnTzK71JOip/HFAqAkWb2mpltNLPP9tXQzJab2QYz21BRUTGIJYqInJo1eys5N39U1BxgPs7PUIgDzgH+Crgc+FczK+itoXPuPudcoXOuMCcnuvrnRGT4KatvZXdZA0umR1fXEfg7JLUEqHTONQFNZrYaOBPY42NNIiKnbc3eSgDOj8JQ8HNP4U/A+WYWZ2YjgIXATh/rEREZEGuKKslKSWDWmHS/Szlpnu0pmNmjwIVAtpmVAHcD8QDOuXudczvN7HlgC9AF/No51+fwVRGRaOCc4429lSyZnk1MjPldzknzLBScczeG0eYHwA+8qkFEZLDtKm2gsrEt6oaiHqczmkVEBtDqPcERkl5f98ArCgURkQG0akcZc8alMyYjye9STolCQURkgJTXt7LxUA3L5ozxu5RTplAQERkgL+woA2DZXIWCiMiwt2p7KVOyU5iWm+p3KadMoSAiMgBqm9t5e18Vl88dg1n0DUU9TqEgIjIAXt5ZTmeXi+rjCaBQEBEZEM9vL2VsRhLzxmf4XcppUSiIiJymprZOVu+p4PI50d11BAoFEZHT9vqeCto6u7g8yruOQKEgInLaVm45yqiUBBbkj/S7lNOmUBAROQ01Te28tKOcq+ePIy42+j9So38NRER8tOK9o7QHurjunAl+lzIgFAoiIqfhjxsPM2dcOrPHRd+1E3qjUBAROUU7j9Wz7Ug9154z3u9SBoxCQUTkFD2+sYT4WOPq+Xl+lzJgFAoiIqegI9DF05uOcPHM0YxKSfC7nAGjUBAROQWv7iqnqqmd6wqHTtcRKBRERE7Jo+uKyU5N5IKC6LzCWl8UCiIiJ6movJFXd1fwmfMmEj8Ezk3oaWitjYjIIPjNmgMkxMVw83mT/C5lwCkUREROQlVjG0++W8Inz84jKzXR73IGnEJBROQkPLy2mLbOLm5bMtnvUjzhWSiY2f1mVm5m2/ppt8DMAmZ2rVe1iIgMhNaOAA++fZCLZuQwLTfN73I84eWewgPAshM1MLNY4HvACx7WISIyIFZsPkplYzu3nz/F71I841koOOdWA9X9NPsi8ARQ7lUdIiIDoTPQxb2v72PW2HQWTc3yuxzP+HZMwczygGuAe8Nou9zMNpjZhoqKCu+LExH5gD9uLGF/ZRNfubQg6q+udiJ+Hmj+MfB151ygv4bOufucc4XOucKcnKF1ooiIRL7WjgA/eWkvZ0/M5JJZuX6X46k4H9+7EHgslLjZwJVm1umce9rHmkREPuTBtw9SWt/Kjz81f0jvJYCPoeCc6x7PZWYPACsVCCISaepaOvjZq/tYWpDDeVOG7rGE4zwLBTN7FLgQyDazEuBuIB7AOdfvcQQRkUjwq9X7qWvp4J8vn+F3KYPCs1Bwzt14Em0/51UdIiKn6khtC79es5+PzRvL3LwMv8sZFDqjWUSkD//vuZ04B9+4YqbfpQwahYKISC/WHahm5ZZj/O3SqYwfOcLvcgaNQkFE5AMCXY7vPLOdsRlJ3LF0qt/lDCqFgojIB/xhw2G2H63nm1fOIjkh1u9yBpVCQUSkB+ccP32liHMmjeTj88b6Xc6gUyiIiPRwoLKJI7UtXHNW3pA/Ua03CgURkR7e3FcFwOJp2T5X4g+FgohID28VVTIuI4n8rOEz4qgnhYKISEhXl+Pt/VUsmpY9LLuOQKEgItJtx7F6aps7WDxt6M9x1BeFgohIyJtFlQAsmjo8jyeAQkFEpNub+6qYlpvK6PQkv0vxjUJBRARo7+xi/YFqFg/hS22GQ6EgIgJsKq6hpSPAomE6FPU4hYKICMGuoxhjWFxI50QUCiIiBA8yn5GXQUZyvN+l+EqhICLDXl1zB5sP17Jk+vDuOgKFgogIr++tINDl+OjM0X6X4juFgogMe6/sLGNUSgLzJ2T6XYrvFAoiMqx1Brp4bU8FF87IITZmeE5t0ZNCQUSGtXeLa6lt7uBidR0BCgURGeZe3lVGXIxxQYEOMoNCQUSGuVd2lrNwyijSkob3UNTjPAsFM7vfzMrNbFsfz99kZltCt7fM7EyvahER6U1xVTN7yxs16qgHL/cUHgCWneD5A8BS59w84B7gPg9rERH5kJd3lQFwyaxcnyuJHHFeLdg5t9rM8k/w/Fs97r4DjPeqFhGR3ryyq5ypOSlMykrxu5SIESnHFG4D/ux3ESIyfNS1dPDO/iounqWuo54821MIl5ldRDAUlpygzXJgOcDEiRMHqTIRGcpe3FFGR8BxxdwxfpcSUXzdUzCzecCvgaudc1V9tXPO3eecK3TOFebk5AxegSIyZK3ccpS8zGSdxfwBvoWCmU0EngRuds7t8asOERl+apraWbO3ko/NG4uZzmLuybPuIzN7FLgQyDazEuBuIB7AOXcv8C0gC/h5aKN0OucKvapHROS4VTtK6exyfGzeOL9LiThejj66sZ/nbwdu9+r9RUT6snLLMSZljWBuXrrfpUScSBl9JCIyKKoa23hrX5W6jvqgUBCRYeXP20oJqOuoTwoFERlWnt1yjCk5Kcwck+Z3KRFJoSAiw0ZZfStrD1TxsXnj1HXUB4WCiAwbj28socvBJ87K87uUiKVQEJFhwTnHHzccZuHkUeRna66jvigURGRYWHegmoNVzVxfOMHvUiKaQkFEhoU/bCghNTGOK88Y63cpEU2hICJDXkNrB89tPcbHzxxHckKs3+VENIWCiAx5K7cco6UjwA0L1HXUH4WCiAx5v19/mILRqZw5PsPvUiJe2KFgZkvM7NbQzzlmNtm7skREBsau0no2H67l+sIJOjchDGGFgpndDXwd+GbooXjgIa+KEhEZKL9dc5Ck+BiuPUdX/A1HuHsK1wBXAU0AzrmjgM4RF5GIVtXYxlObj/DJs8eTOSLB73KiQrih0O6cc4ADMDOd+SEiEe+RtcW0d3Zx6+J8v0uJGuGGwh/M7JdAppn9DfAS8CvvyhIROT3tnV387zuHuKAgh2m56tgIV1gX2XHO/aeZXQrUAzOAbznnXvS0MhGR0/Dc1mOUN7Tx/Wvz/S4lqoQVCqHuoleccy+a2QxghpnFO+c6vC1PROTkOee4/80DTM1J4YLpOX6XE1XC7T5aDSSaWR7BrqNbgQe8KkpE5HSsO1DNlpI6Prd4MjExGoZ6MsINBXPONQOfAP7HOXcNMNu7skRETt1PXy0iOzWBa8/WMNSTFXYomNlHgJuAZ0OPhdX1JCIymDYfruWNvZXcfv4UzXN0CsINhX8AvgE86ZzbHjqb+RXvyhIROTU/fWUvmSPi+cx5k/wuJSqF+22/GegCbjSzzwBG6JwFEZFIsf1oHS/tLOcrlxaQmqjOjFMR7m/tYeCfgG0Ew6FfZnY/8DGg3Dk3t5fnDfgJcCXB0Pmcc+7dMOsREfmQn71aRFpiHLcsyve7lKgVbvdRhXPuGefcAefcoeO3fl7zALDsBM9fAUwP3ZYDvwizFhGRD9l5rJ4/byvllkX5ZCTH+11O1Ap3T+FuM/s18DLQdvxB59yTfb3AObfazPJPsMyrgQdD02e8Y2aZZjbWOXcszJpERABo6wzwlT+8x6gRCdy2RBM4n45wQ+FWYCbB2VGPdx85oM9QCEMecLjH/ZLQYwoFETkpP1y1h53H6vnNLYWMTNHEd6cj3FA40zl3xgC/d29nlPR68NrMlhPsYmLixIkDXIaIRLM3iyq5b/V+blo4kYtnjfa7nKgX7jGFd8xsoE9WKwF6XhtvPHC0t4bOufucc4XOucKcHJ2yLiJBtc3tfPUP7zElJ4V/+SudTzsQwg2FJcBmM9ttZlvMbKuZbTnN914BfNaCzgPqdDxBRMLVEejizkc2UdXUxk9uOEsnqg2QcLuPTjSKqFdm9ihwIZBtZiXA3QSPSeCcuxd4juBw1CKCQ1JvPdn3EJHh656VO1hTVMn3r53HGbr28oAJd+rs/oaf9vaaG/t53gF/f7LLFRF58O2DPPj2IZZfMIXrCyf0217CF273kYhIRHhtdznfeWYHF8/M5evLZvpdzpCjUBCRqLF2fxVfeGgjBaPT+MmNZxGrabEHnEJBRKLCpuIaPv/AesaPHMH/3nau5jbyiEJBRCLetiN13HL/OrLTEnn49oVkpyb6XdKQpVAQkYj29r4qbrzvHdKS4nn49oWMTk/yu6QhTaEgIhHr+W3HuOW36xiTkcQfv/ARxo8c4XdJQ5465UQk4jjn+N1bB/nuyh3Mn5DJ/Z9bQOYIzWk0GBQKIhJRWtoD3PXUVp7cdIRLZo3mf27U2cqDSaEgIhGjuKqZLzy0kZ2l9Xz5kgK++NFpxGjY6aBSKIiI75xzPLXpCN/603ZiDO6/ZQEXzcz1u6xhSaEgIr6qbW7nrqe38eyWYyzIH8mPrp/PhFE6oOwXhYKI+GZvWQO3PrCe0rpWvrZsBn97wVSdpewzhYKI+GL1ngr+/uF3SUqI5fE7FjF/QqbfJQkKBREZZM45HnrnEN9+ZgcFo9P4zS2FjMtM9rssCVEoiMigaesMcPeftvPY+sNcPDOXn9x4luYwijDaGiIyKMrqW/nCQxvZVFzLnRdN48uXFuj4QQRSKIiI53aV1nPL/etoaO3kFzedzRVnjPW7JOmDQkFEPLXuQDW3/W49IxJieeKORcwam+53SXICCgUR8cyq7aXc+egmxo9M5sHPn6sJ7aKAQkFEBpxzjt++eZB/e3YHZ4zP5LefW8CoFE1oFw0UCiIyoDoDXXz7me089E4xl80ezY8/NZ8RCfqoiRbaUiIyYKqb2vmHxzbxxt5K/nbpFL5++UxNaBdlFAoiMiA2HKzmi49uoqqxne998gxuWDDR75LkFCgUROS0BLocv1y9jx+u2sP4kck8+XeLmJuX4XdZcoo8vRynmS0zs91mVmRm3+jl+Ylm9qqZbTKzLWZ2pZf1iMjA2lfRyPW/fJvvP7+bZXPGsPKLSxQIUc6zPQUziwV+BlwKlADrzWyFc25Hj2b/AvzBOfcLM5sNPAfke1WTiAyMzkAXv1lzgB+9uIek+Fj+64Yz+ev5eZjp+EG087L76FygyDm3H8DMHgOuBnqGggOOn8mSARz1sB4RGQDv7K/i7j9tZ3dZA5fNHs2/XTOX3LQkv8uSAeJlKOQBh3vcLwEWfqDNt4FVZvZFIAW4pLcFmdlyYDnAxIk6eCXih6O1LXzv+V38afNR8jKTufcz53D5nNHaOxhivAyF3v6nuA/cvxF4wDn3QzP7CPC/ZjbXOdf1vhc5dx9wH0BhYeEHlyEiHqpr6eDnrxXx2zcPAvCli6dzx9KpJCfE+luYeMLLUCgBJvS4P54Pdw/dBiwDcM69bWZJQDZQ7mFdIhKGQ1VNPLb+MI+uK6aupYNrzsrjq5fNIE/XPhjSvAyF9cB0M5sMHAE+BXz6A22KgYuBB8xsFpAEVHhYk4icQEegi5d2lPHw2mLWFFUSG2NcPDOXf7ykgNnjNJHdcOBZKDjnOs3sTuAFIBa43zm33cy+C2xwzq0Avgr8ysy+TLBr6XPOOXUPiQyy6qZ2Hnz7II+uK6asvo28zGS+cmkB1xdOYEyGDiIPJxZtn8GFhYVuw4YNfpch4ql3i2v4yUt7yRuZzJJp2XxkShYjPZhQrjPQxSPrivnPF3ZT39rJ0oIcbj5vEhfNzNUFcIYYM9vonCvsr53OaBaJIIEux72v7+NHL+5h5IgENh6q4ZG1xZjB7LHpLJqaxaKp2Zw7eRQpp3kZy7f2VfLdZ3awq7SBRVOzuPvjc5gxJm2A1kSilUJBJEI0tHaw/MGNvL2/io+fOY5/v2YuI+Jjea+kjreKKnlzXyW/e+sQv3rjAHExxlkTM1k0NZulM3I4c3xm2N/stx+t43vP72b1ngrGZSTx85vO5oq5YzS0VAB1H4lEjH95eiuPrC3mPz45j+vOGd/rh3RrR4ANB2t4c18lbxVVsvVIHV0OslISuGhmLpfPGcPSghwS4j48g83h6mZ+uGo3T28+SkZyPH9/0VQ++5F8kuI1tHQ4UPeRSBR5t7iGh9cWc+uiyVxfOKHPdknxsSyZns2S6dkA1DV38PreCl7eWcaLO8p4fGMJ6UlxXHnGWBbkjyJzRDzpyfG8sK2UB98+hBncceFUvrB0KhnJ8YO1ehJFtKcg4rOOQBcf/5811LV08OJXlpJ6iscKOgJdrCmqZMXmo6zaXkpTe6D7OTO47pzxfPnSAsZm6DyD4Uh7CiJR4v41B9hV2sAvbz7nlAMBID42hotm5HLRjFxaOwKU1bdS29xBTXM7E0eNYEpO6gBWLUOVQkHER4erm/nxS3u5dPZoLp8zZsCWmxQfy6SsFCZlDdgiZZjw9HoKItI35xz/56mtxBh8+6o5fpcjAigURHzzxLtHeGNvJV+/YqbmE5KIoVAQ8UFFQxv3rNxB4aSRfGbhJL/LEemmUBDxwbdXbKelPcB/fHIeMZpOQiKIQkFkkD2/7RjPbj3GP1wynWm5GhEkkUWhIDKIyutb+eaTW5mbl87yC6b4XY7IhygURAaJc45/enwLLR0BfnzDWcTH6s9PIo/+V4oMkgffPsTqPRXc9Vez1W0kEUuhIDII9pY18H+f28lFM3L4zMKJfpcj0ieFgojHWjsC3PnIJlIT4/jetfM0RbVENE1zIeKx7zyzg91lDfzu8+eSm6ZLW0pk056CiIeeee8oj64r5o4Lp7K0IMfvckT6pVAQ8cihqia++eRWzpk0kq9cWuB3OSJhUSiIeKCxrZPlD24kNsb47xs1/FSih44piAywQJfjHx/bRFFFIw9+/lxNdidRRV9fRAbYf67azUs7y7n747NZPC3b73JEToqnoWBmy8xst5kVmdk3+mhzvZntMLPtZvaIl/WIeO2JjSX84rV9fHrhRG4+T7OfSvTxrPvIzGKBnwGXAiXAejNb4Zzb0aPNdOCbwGLnXI2Z5XpVj4jXXtlVxtee2MKiqVl856o5Oh9BopKXewrnAkXOuf3OuXbgMeDqD7T5G+BnzrkaAOdcuYf1iHhm3YFq7njoXWaPTeeXN5+jA8sStbz8n5sHHO5xvyT0WE8FQIGZvWlm75jZst4WZGbLzWyDmW2oqKjwqFyRU7PhYDW3PbCevJHJPHDrAtKS4v0uSeSUeTn6qLd9Z9fL+08HLgTGA2+Y2VznXO37XuTcfcB9AIWFhR9chsiga2rr5Nmtx3hkbTGbD9cyLiOJh25bSFZqot+liZwWL0OhBJjQ4/544Ggvbd5xznUAB8xsN8GQWO9hXSInraKhjfUHq1l3oJqNh2rYcayeQJdjak4K//qx2Vx79ngyRmgPQaKfl6GwHphuZpOBI8CngE9/oM3TwI3AA2aWTbA7ab+HNYn0yzlHSU0L6w9Ws/5gNWsPVLO/ogmApPgY5k/I5I6lU7mgIIcF+SN1QFmGFM9CwTnXaWZ3Ai8AscD9zrntZvZdYINzbkXoucvMbAcQAP7ZOVflVU0ivaluamfXsXq2HKljc3Etmw/XUlrfCkBaUhwL8kdxQ+EEFk7JYs64dB1EliHNnIuuLvrCwkK3YcMGv8uQKNTaEWBvWSO7SuvZXdrA7rIGdpU2UNHQ1t1mUtYI5o3PpHDSSBbkj2LGmDRiY7QnINHPzDY65wr7a6dpLmRIamrrZMexeraW1LHtaB07jtZTVN5IZ1fwS1BiXAzTR6dywfQcZo5JY8aYNM7Iy2BkSoLPlYv4S6EgUa+2uZ31B2vYVFzDnrJG9pY3UFzdzPGd4Jy0ROaMS+fiWbnMHpvBzLFp5GelaA9ApBcKBYk6rR0B1h6o5rXd5bxVVMXusgYA4mKMydkpzB2XwSfOGs/cvHTOyMsgN10XthEJl0JBokJ1Uzsv7Sxj1fYy1hRV0NrRRUJcDAsnj+Kq+eNYkD+KeeMzSIqP9btUkaimUJCI1Rno4uVd5Tyytpg39lbQ5WBcRhLXF07gohm5nDcli+QEhYDIQFIoSMTZXdrAiveO8PjGEsrq2xiTnsQdF05l2ZyxzM1L13kBIh5SKIjvnHPsLW9k1fZSVm45xq7SBmIMzp+ewz1XT+SjM3OJ07kBIoNCoSC+qG/tYP2Bat7aV8Uru8o5UBk8Y/jsiZl856o5XHnGWHLSNI+QyGBTKIjnjk8b8W5xDe8equHd4lq2H62jy9F9sPi2JZO5dPZoRmukkIivFAoy4Gqb29lV2sDOY/VsPFTD+oPVlNUHzxoekRDLvPEZ3HnRNM6bmsXZE0dqxJBIBFEoyElr7QhQ0dDG0doWjta1cLS2lUNVTRysauZAZdP7po0Ym5HEwslZLMgfyTmTRlEwOlXHB0QimEJBcM5R39JJWUMrpXWtVDS0UdPcTk1zO9VNHdQ0tVPd1E5VUxvlDW00tHZ+aBk5aYnkZ41gaUEOBaNTmTEmnZlj0tQdJBJlFArDQGtHgKLyRnaVNlBS00xVY/ADvqKhjbL6NsobWmnt6PrQ62JjjMzkeEalJDAqJYGC0WksmZZNTloiuWlJjM1MYmxGMuMykxiRoP9KIkOB/pKHkM5AF4drWigqb2RPWbBPf3dpA/srmwh0/WU23MwR8WSlJJCdmshZEzPJTUtkdHpS9y0nLZFRIxJIS4ojRvMDiQwrCoUocrybp6S2mSM1LRypbeFQVTPF1c0cqmqiuLqZjsBfPvzHj0xm5ph0Lp8zhllj05k5No2Jo0boegAi0ieFQoRxznG0rpW9ZQ0cqGzqvh2tbeFYXSvN7YH3tR+REMvEUSOYlpvKZXPGMCU7hSk5qUwfnUq6LiAvIidJoeAT5xzF1c3sPNbAvopG9lc0sa+ikaLyRhrb/nIgNy0xjsk5KRSMTmNpQS7jMpPIy0wmb2Qy4zKTyUpJ0LQPIjJgFAoec85R3tDG/oomisob2FveyO7SBnYcraehx4d/bloiU3JS+MTZeUwfnUZBbipTc1P1oS8ig0qhcBqcczS2db5vFM+xutbg+P3aFg5Xt3Couul9I3tSE+OYPjqVq88ax9xxGcwam86UnBTS1NUjIhFgWIeCc472QBct7QGa2gM0t3XSGLo1tHbS2NpJXUtH9622pYPa4+P3G9upbGqnvfPDQznTk+IYl5nMhFHJLJmeTX7WCPKzU5iWm8qY9CR98xeRiDVsQuHV3eXcs3IHbR1dtHUGaO3ooqUj8L6hmn2JMUhPjiczOZ7MEQnkpCYyY3Q62akJZKUGh3YGh3MG/9W3fhGJVsMmFDKS45k1Jp3E+BiS4mNJiotlREIsyQmxJMfHkpIYS0piHCkJcaQmxZGWFEdqYhwZyfGkJGi8vogMD8MmFM6eOJKzbxrpdxkiIhHN07OYzGyZme02syIz+8YJ2l1rZs7MCr2sR0RETsyzUDCzWOBnwBXAbOBGM5vdS7s04EvAWq9qERGR8Hi5p3AuUOSc2++cawceA67upd09wPeBVg9rERGRMHgZCnnA4R73S0KPdTOzs4AJzrmVHtYhIiJh8jIUehuu0z3+08xigP8CvtrvgsyWm9kGM9tQUVExgCWKiEhPXoZCCTChx/3xwNEe99OAucBrZnYQOA9Y0dvBZufcfc65QudcYU5Ojocli4gMb16GwnpguplNNrME4FPAiuNPOufqnHPZzrl851w+8A5wlXNug4c1iYjICXgWCs65TuBO4AVgJ/AH59x2M/uumV3l1fuKiMipM+f6n+YhkphZBXDI47fJBio9fo/BonWJPENlPUDrEql6W5dJzrl++9+jLhQGg5ltcM4NiRPptC6RZ6isB2hdItXprIuuyygiIt0UCiIi0k2h0Lv7/C5gAGldIs9QWQ/QukSqU14XHVMQEZFu2lMQEZFuCgXAzO4xsy1mttnMVpnZuD7a3WJme0O3Wwa7znCY2Q/MbFdofZ4ys8w+2h00s62hdY7IEwZPYl3CmqLdL2Z2nZltN7OuE00PHyXbJNx1iehtAmBmo8zsxdDf84tm1usFV8wsENomm81sRW9t/NDf79jMEs3s96Hn15pZflgLds4N+xuQ3uPnLwH39tJmFLA/9O/I0M8j/a69lzovA+JCP38P+F4f7Q4C2X7Xe7rrAsQC+4ApQALwHjDb79o/UOMsYAbwGlB4gnbRsE36XZdo2CahOr8PfCP08zdmg0MKAAAFhUlEQVRO8LfS6Hetp/I7Bv7u+GcZwRklfh/OsrWnADjn6nvcTaHHxH09XA686Jyrds7VAC8CywajvpPhnFvlgmeTQ3DqkPF+1nM6wlyXcKdo941zbqdzbrffdQyEMNcl4rdJyNXA70I//w74ax9rOVnh/I57rt/jwMVm1u91hRUKIWb272Z2GLgJ+FYvTfqdCjwCfR74cx/POWCVmW00s+WDWNOp6mtdonG79CXatklfomWbjHbOHQMI/ZvbR7uk0CzN75hZpARHOL/j7jahL1d1QFZ/Cx4212g2s5eAMb08dZdz7k/OubuAu8zsmwTnbLr7g4vo5bW+DN3qb11Cbe4COoGH+1jMYufcUTPLBV40s13OudXeVNy3AViXiNgu4axHGKJmm/S3iF4ei7i/lZNYzMTQdpkCvGJmW51z+wamwlMWzu/4lLbDsAkF59wlYTZ9BHiWD4dCCXBhj/vjCfarDrr+1iV0EPxjwMUu1KHYyzKOhv4tN7OnCO6ODvoH0ACsS39TtA+Kk/j/daJlRMU2CUNEbBM48bqYWZmZjXXOHTOzsUB5H8s4vl32m9lrwFkE+/P9FM7v+HibEjOLAzKA6v4WrO4jwMym97h7FbCrl2YvAJeZ2cjQKIXLQo9FFDNbBnyd4DTkzX20SQldGxszSyG4LtsGr8rwhLMu9DNFe7SIlm0SpmjZJiuA46MIbwE+tBcU+ntPDP2cDSwGdgxahX0L53fcc/2uBV7p60vi+/h9FD0SbsATBP8AtwDPAHmhxwuBX/do93mgKHS71e+6+1iXIoL9iJtDt+OjD8YBz4V+nkJwtMJ7wHaC3QK+134q6xK6fyWwh+C3t4hbF+Aagt/a2oAy4IUo3ib9rks0bJNQjVnAy8De0L+jQo93/90Di4Ctoe2yFbjN77pP9DsGvkvwSxRAEvDH0N/ROmBKOMvVGc0iItJN3UciItJNoSAiIt0UCiIi0k2hICIi3RQKIiLSTaEgw4aZNZ7m6x8PndV6ojavnWj20HDbfKB9jpk9H257kdOhUBAJg5nNAWKdc/sH+72dcxXAMTNbPNjvLcOPQkGGHQv6gZltC12/4IbQ4zFm9vPQ9QJWmtlzZnZt6GU30eOMVzP7RWiStO1m9p0+3qfRzH5oZu+a2ctmltPj6evMbJ2Z7TGz80Pt883sjVD7d81sUY/2T4dqEPGUQkGGo08A84EzgUuAH4TmvvkEkA+cAdwOfKTHaxYDG3vcv8s5VwjMA5aa2bxe3icFeNc5dzbwOu+fTyvOOXcu8I89Hi8HLg21vwH47x7tNwDnn/yqipycYTMhnkgPS4BHnXMBoMzMXgcWhB7/o3OuCyg1s1d7vGYsUNHj/vWh6a3jQs/NJjhNSk9dwO9DPz8EPNnjueM/byQYRADxwE/NbD4QAAp6tC8nOJWEiKcUCjIc9XWhkRNdgKSF4FwymNlk4J+ABc65GjN74Phz/eg5p0xb6N8Af/k7/DLB+YTOJLgX39qjfVKoBhFPqftIhqPVwA1mFhvq57+A4IRha4BPho4tjOb9U6XvBKaFfk4HmoC6ULsr+nifGIKzUwJ8OrT8E8kAjoX2VG4meMnF4wqI3llTJYpoT0GGo6cIHi94j+C3968550rN7AngYoIfvnuAtQSvVgXBa2xcCLzknHvPzDYRnM10P/BmH+/TBMwxs42h5dzQT10/B54ws+uAV0OvP+6iUA0intIsqSI9mFmqc67RzLII7j0sDgVGMsEP6sWhYxHhLKvROZc6QHWtBq52weuDi3hGewoi77fSzDKBBOAe51wpgHOuxczuJnjd2+LBLCjUxfUjBYIMBu0piIhINx1oFhGRbgoFERHpplAQEZFuCgUREemmUBARkW4KBRER6fb/Ae2uq9iNqE4KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2121962d550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lasso picked 25 features and eliminated the other 7 features\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEICAYAAAD7pTujAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xu8VVW99/HP1ysqJqlkkiIJGikoBmqmFF7SjpqXEz54ycRMM/X4dDpeeLSMSk8WlnlJLQ2vZKapR8VEwhQkRECBLYF3PZKkaHnLJJXf88cYWxaLfVlr7bXXWnvzfb9e+8VaY44552+uvdljjzHH/A1FBGZmZo1kjXoHYGZmVsyNk5mZNRw3TmZm1nDcOJmZWcNx42RmZg3HjZOZmTUcN07WpUn6hKRHJb0p6VRJ60m6U9Lrkm6WdJSke0s4zlmSrqpFzG3E0FfSW5LWrNLxxkq6oRrHWt1Iek7SPiXU6ycpJK1Vi7hWJ26crCYkHSlpdv7lu0TS7yXtUYVDnwHcHxEbRsTFwEhgM2CTiDgsIiZExL7tHSQi/jsivtbRYDryyyoi/jciekbE+xWcd4SkxeXuV8bxr5F0bmcdv1I5rpB0UFH5z3L56DqFZh3kxsk6naRvAT8D/pvUcPQFLgMOrsLhtwIWFL1/IiLeq8KxrWt4Ajim+U3+w+Aw4Om6RWQd5sbJOpWkjYDvAydHxK0R8Y+IeDci7oyI03OddfNfui/mr59JWrfgGAdKmivpNUl/krRDLr8P2BO4NPfIbgTOAUbl98dJGi3pwYJjbS9psqS/SXpJ0lm5fKUhMEmfzud6TdI8SSMKtt0v6QeSpufhxHslbZo3T83/vpZj2E3SAEkP5KHGVyTd1MpntVKvq53zFO63AfB7oE8+51uS+uTN60i6Lu+/QNKwgv36SPqdpKWSnpV0ainf0xbOf5GkFyS9IWmOpOEF23bJPeY38uf901zeQ9INkl7Nn/EsSZsVxHVH/h49Jen4dkK4E9hd0ofz+y8A84G/FsSxhqRvS3pe0sv5M9moYPvRedurks4uur41JI2R9HTe/ltJG1fyWVnp3DhZZ9sN6AHc1kads4FPA0OAHYFdgG8DSPoUMB74OrAJ8AvgDknrRsRewDTglDwcdgSpd3ZTfv+rwpNI2hD4A3AP0AcYAEwpDkbSx4CJwLnAxsBpwO8k9S6odiRwLPARYJ1cB+Cz+d9eOYYZwA+Ae4EPA1sAl7TxWRRr7TwfiIh/AP8GvJjP2TMiXsybDwJ+A/QC7gAuzde4BumX+jzgY8DewDcl7VdGbM1mkb53GwO/Bm6W1CNvuwi4KCI+BPQHfpvLjwE2ArYkfV9PBP6Zt90ILCZ9j0YC/y1p7zbO/06+tsPz+68A1xXVGZ2/9gS2Bnqy4rPYDrgcODqfcxPS96nZqcAhwOfy9r8DP28jHqsCN07W2TYBXmlnmO0o4PsR8XJELAW+R/pFAXA88IuImBkR70fEtcAyUmNWrgOBv0bETyLinYh4MyJmtlDvy8DdEXF3RCyPiMnAbGD/gjpXR8QTEfFP0i/cIW2c913ScGOffN4H26hbrJzztOTBfB3vA9eTGn+AnYHeEfH9iPhXRDwDXMmKX/Ali4gbIuLViHgvIn4CrAt8Im9+FxggadOIeCsiHioo3wQYkL+vcyLiDUlbAnsAZ+bPai5wFSt+HlpzHfCV3Bv6HHB70fajgJ9GxDMR8Rbw/4DDcy91JHBXREyNiGXAd4DlBft+HTg7Ihbn7WOBkfIkiE7lxsk626vApu38R+4DPF/w/vlcBumX+n/loZ/XJL1G+mu7D+XbktLuQ2wFHFZ0zj2AzQvq/LXg9dukv8RbcwYg4OE8tPbVMmIu5zyl7N8jfy+2Ig0DFl7jWaR7gmWR9F+SFuZhy9dIPaLm4cfjgG2BRXno7sBcfj0wCfiN0lDujyWtTfq+/i0i3iw4xfOk3l2rcoPfm9Tjvis35oVa+hlbK19vH+CFgmP9g/Rz22wr4LaCz2kh8D4VfFZWOrf81tlmkIZdDgFuaaXOi6w8saFvLoP0S+O8iDivCrG8ABxRYr3rI6K9ex0tWSXNf0T8ldQDRGmG4h8kTY2Ipyo4fsnnbccLwLMRsU1HTprvL51JGhZcEBHLJf2d1BgTEU8CR+RhxH8HbpG0SW4Avgd8T1I/4G7gcdLw58aSNixooPoCfykhnBtI9xz3bGFb889Ys77Ae8BLwBLgkwXXtD6pV9fsBeCrETG9hevvV0JcVgH3nKxTRcTrpF8YP5d0iKT1Ja0t6d8k/ThXuxH4tqTe+Yb/OaRfNJCGmk6UtKuSDSQdkO8flesu4KOSvqk0CWNDSbu2UO8G4IuS9pO0Zr55P0LSFi3ULbaUNCS0dXOBpMMK9v07qSEpe7p4O14CNim8yd+Oh4E3JJ2p9GzYmpIGSdq5jX2aP4vmr3WADUm/5JcCa0k6B/hQ8w6Sviypd0QsB17Lxe9L2lPSYKVnut4gDfO9HxEvAH8CfpjPsQOp9zWhhGu6GPg8KyalFLoR+E9JH5fUkxX3Jt8j/dF0oKQ98jV9n5V/N14BnCdpq3xNvSVVY6aptcGNk3W6iPgp8C3SkMtS0l+ip7DivsC5pHs684Em4JFcRkTMJvU6LiX9Yn+KdGO7kjjeJP3y+iJpuOtJWvgrO/+CPJg0zNUc7+mU8P8lIt4GzgOm52GgT5Pu78yU9Bbpxv3/jYhnK7mGNs67iPQL+Jl83jaHPfM9qC+S7mE9C7xCurfTVuM2hjRpofnrPtLQ3O9J07mfJ/WSXyjY5wvAgnztFwGHR8Q7wEdJjcIbpGGyB1jxB8kRQD9Sb+c24Lv5vl97n8HfImJKtLxI3XjSUOLUfL3vAP+R91sAnEyazLGE9HNW+MzYRaTv272S3gQeAlr6o8aqSF5s0MzMGo17TmZm1nDcOJmZWcNx42RmZg3HjZOZmTUcP+dUoU033TT69etX7zDMzLqUOXPmvBIRvdur58apQv369WP27Nn1DsPMrEuR9Hz7tTysZ2ZmDcg9JzOriin39a93CFYje+/V+UtldZuek6Rekk4qeD9C0l1l7H+N0po2c/NXudmfzcysSrpN40Rar+akdmu17fSIGJK/5lYjKDMzK19DNU5KK4EuknSVpMckTZC0j9JKoE8qrao5VtJ4pVVCn9GK1TvPB/rnXs+4XNZT0i35mBMkqYPxnaC0qufspUuXduhazcysdQ3VOGUDSIkWdwAGklYC3YO0AuhZuc5AYD/SiqnfzevAjAGezr2e03O9nYBvAtuRskTv3s65z5M0X9KFKlgmvFlE/DIihkXEsN69250JaWZmFWrExunZiGjKKfYXAM1ZhptImYoBJkbEsoh4BXiZ1hf9ejivXrkcmFuwf0v+H6nR25m03PSZHb4SMzOrSCPO1ltW8Hp5wfvlrIi3sM77tH4dpdYjIpY07yPpalJPzcxKVIsZXLb6aMSeU6XeJC18VhFJm+d/RVq19bEqxWVmZmVqxJ5TRSLi1Txx4jHS4mcTyzzEBEm9SctLzwVOrHaMZmZWGi82WKFhw4aF0xeZmZVH0pyIGNZeve40rGdmZt1Elx3Wk9QLODIiLitzv9uAjxcVnxkRk6oWnNlqaOzYsfUOoVtZ3T/PLts4sSIjRFmNU0Qc2jnhmJlZtXTlYb2VMkJIOl3SrPwQ7fegtIwTud5YSddLui+XH1/XKzMzW8115cbpg4wQwGRgG1LGiCHAUEmfzfVKyThB3n4AsBtwjqQ+xSd0+iIzs9royo1ToX3z16PAI6RGaJu8rZSMEwD/ExH/zFkn/khq6Fbi9EVmZrXRle85FRLww4j4xUqFUj9KyzgBUDyn3nPszcqwut/At+rqyj2nwowQk4CvSuoJIOljkj5S5vEOltRD0ibACGBW1SI1M7OydNmeUwsZIX4NzMirYrwFfJmUT69UD5OySvQFfhARL1Y5ZDMzK1GXbZwAIuLIoqKLWqg2qKD+6ILXzxVuA56IiBOqGZ+ZmVWmKw/rmZlZN9Wle07VEhFj6x2DmZmt0G0ap+J0RpJGAKdFxIEl7i/gXOAw0r2qyyPi4k4K16zbWTxmWr1D6PK2OH94vUNoGN1pWK85nVGlRgNbAgMj4pPAb6oRlJmZla+hGqdS0g3lVEPjJd0v6RlJp+bdV0pnlMt6SrolH3NC7h215hvA9/PDukTEy514qWZm1oaGapyyUtINDQT2I2Vx+K6ktSlIZxQRp+d6OwHfBLYDtgZ2b+O8/YFROT3R7yVtU1zB6YvMzGqjERunUtINTYyIZTnV0MvAZq0c6+GIWJyPNZeV0xUVWxd4Jy+CdSUwvriC0xeZmdVGI06IKCXdUGGd92n9OkqtB7AY+F1+fRtwdSnBmlnim/lWTY3Yc6pUYTqjStwO7JVffw54osMRmZlZRRqx51SRFtIZTSzzEOcDEyT9Jyn90deqHaOZmZVG6XaOlWvYsGExe/bseodhZtalSJqT7+23qTsN65mZWTfRbYb1SiXpNuDjRcVnRsSkesRjZmar6jaNU6npiyLi0Fb23wu4AFgHmAMcFxHvdWrQZt3IT0aVlCmsW/qvm+6qdwjdTnca1qs4fZGkNYBrgcMjYhDwPHBMFWMzM7MyNFTjVMf0RZsAyyKiefr4ZOBLLcTnDBFmZjXQUI1TVo/0Ra8Aa0tqnkEykpQEdiXOEGFmVhuN2DjVPH1RPv7hwIWSHiY90Ov7TWZmddKIEyLqkr4oImYAwwEk7QtsW3rIZuZJAVZNjdhzqlSH0hdJ+kj+d13gTOCKKsVlZmZl6jaNU0S8CkzPEynGtbvDqk6XtBCYD9wZEfdVN0IzMyuV0xdVyOmLzMzK5/RFZmbWZTXihIhO5fRFZmaNr0s1TqWmKGpj/1OAwaQl2XvnqegouRjYH3gbGB0Rj3TCJZh1Wz8/sfvcpj35ir3ar2SdqqsN61WcoiibDuxDSk9U6N+AbfLXCcDlHTiHmZl1UM0bpzqmKCIiHo2I51rYdDBwXSQPAb0kbd5C7E5fZGZWA/XqOdUjRVFbPga8UPB+cS5bidMXmZnVRr0ap5qnKGpHS70tz7E3M6uTek2IqEuKojYsZuVEr1sAL1ZwHLPVlicRWDV1tQkRHUpR1IY7gK/kWXufBl6PiCWdcB4zMytBl2qcOpqiSNKpkhaTekbzJV2VN90NPAM8BVxJx2YEmplZBzl9UYWcvsjMrHxOX2RmZl1Wl8oQUSqnKDKrvYUDP1nvEFr1yUUL6x2ClakeD+GOltSn4P1zkjat5jki4lDSbLsR+euX7TVMku6R9Jokr5hmZlZn9RjWGw30aa9SKSS1tbLt/hHxGqWnPBoHHF2NuMzMrGPabZwkndGcPkjShZLuy6/3lnSDpH0lzZD0iKSbJfXM28+RNCvPrPtlnqY9EhgGTMgpiNbLp/mPvH+TpIF5/w1yCqNZkh6VdHAuH53Pcydwr6TNJU3Nx3tMUvNS6809spZSHq0iIqaQpqq39Vk4fZGZWQ2U0nOaCgzPr4eRctmtTUo31AR8G9gnIj4FzAa+leteGhE7R8QgYD3gwIi4Jdc5Kqcg+meu+0re/3JSCiOAs4H7ImJnYE9gnKQN8rbdgGMiYi9S6qNJETEE2JGUJaJQSymPKuL0RWZmtVHKhIg5wFBJG5KyMTxCaqSGkx5e3Y707BHAOsCMvN+eks4A1gc2JqUpurOVc9xacK5/z6/3BQ6S1NxY9QD65teTI+Jv+fUsYHxuMG+PiOLGyczMuph2G6eIeFfSc8CxwJ+A+aSeTH/gWVJDcUThPpJ6AJcBwyLiBUljSY1La5pTEBWmHxLwpYh4vOjYuwL/KIhvqqTPAgcA10saFxHXtXddZlZdnhFn1VTqhIippOG2qcA04ETS8NlDwO6SBgBIWl/StqxoiF7J96BGFhyr1BREk0j3opSPvVNLlSRtBbwcEVcCvwI+VVSls1IemZlZJym1cZoGbA7MiIiXgHeAaRGxlDT77kZJ80mN1cA8S+5K0j2p20lDb82uAa4omhDRkh8Aa5PSDD2W37dkBDBX0qPAl0hLcXyg1JRHkqYBNwN7S1osab82YjMzs07k9EUVcvoiM7PyOX2RmZl1Wd0yfVFrJA0Gri8qXhYRu9YjHrPuZPC1g+sdQquajmmqdwhWptWqcYqIJmBIW3Uk3QFsnZ/PMjOzOvCwXgFJ/w68Ve84zMxWdzVrnHI6oomS5uWZc6MkDZX0gKQ5kiZJ2jzXPT6nLZon6XeS1s/lh+V950mamst6SLo6pz56VNKeuXy0pFtzQtcnJf24nfh6krJbnNtGHacvMjOrgVr2nL4AvBgRO+Yhs3uAS4CRETEUGA+cl+vemlMf7QgsBI7L5ecA++Xyg3LZyQARMRg4Arg2PwQMaQhvFDAYGCVpyzbi+wHwE+Dt1io4fZGZWW3UsnFqAvaR9KOcnHVLYBAwWdJcUo6+LXLdQZKmSWoCjgK2z+XTgWskHQ+smcv2IE9yiIhFwPPAtnnblIh4PSLeAf4MbNVSYJKGAAMi4rbqXa6ZmVWqZhMiIuIJSUOB/YEfApOBBRGxWwvVrwEOiYh5kkaTHrQlIk7M6YsOID14O4SU5qg1ywpeF6ZGKrYbKX/gc7nORyTdHxEjSrs6M/OMOKumWt5z6gO8HRE3ABcAuwK9Je2Wt68tqbmHtCGwJCdzPargGP0jYmZEnAO8Qup9TW2uk1Mn9QVWysfXnoi4PCL6REQ/Uk/sCTdMZmb1U8up5INJy14sB94FvgG8B1wsaaMcy89I2cu/A8wkDdE1sSI33jhJ25B6S1OAecAiUjqkpny80RGxLKfkMzOzLsjpiyrk9EVmZuVz+iIzM+uyVqsMEQCSZgLrFhUfnbNHmFmlxm7Uycd/vXOPbw2l2/ScJPWSdFLB+xGS7iquFxG75iXbC7+aCva7RJKzRJiZ1VG3aZyAXsBJ7dZqg6Rh+ThmZlZHDdU4SeonaZGkq3KaogmS9pE0Pacg2kXSWEnjJd0v6RlJp+bdzwf650UMmxcV7CnplnzMCWpjCp+kNYFxwBlt1HH6IjOzGmioxikbQFrNdgdgIHAk6dmj04Czcp2BwH7ALsB38/NQY4Cn8zDd6bneTsA3ge2ArYHd2zjvKcAdEbGktQpOX2RmVhuN2Dg9GxFNEbGc9MzTlEjz3ZuAfrnOxIhYFhGvAC8Dm7VyrIcjYnE+1tyC/VeSHxA+jJTrz8zM6qwRZ+sVphxaXvB+OSviLTUtUan1diL12J7KI3/rS3oqIgaUEbfZ6s2z6ayKGrFxqtSbrMgkUZaImAh8tPm9pLfcMJmZ1U8jDutVJCJeBabniRTj2t3BzMwaltMXVcjpi8zMyuf0RWZm1mV1p3tOJZF0G/DxouIzI2JSPeIx6y76jZnYYvlz5x9Q40isO2i4nlN+uLbNLp+k0ZIureT4EXFocfoiYHNJS/MDvHMlfa2i4M3MrCpWu55TG26KiFPqHYSZmVWh5yTpjOYUQpIulHRffr23pBsk7StphqRHJN0sqWfePlTSA5LmSJokafOi464h6VpJ5+b3x0p6QtIDFGR6kPRFSTMlPSrpD5I2y/s+Kal3wbGekrRpB6/V6YvMzGqgGsN6U4Hh+fUwUj67tUkph5qAbwP7RMSngNnAt/L2S4CRETEUGA+cV3DMtYAJpOXSv50bru+RGqXPk9IRNXsQ+HRE7AT8BjgjZ4S4gRVLvO8DzMsZJVrzJUnzcy6+LVuq4PRFZma1UY1hvTnAUEkbkjIyPEJqpIYDd5Aakuk588I6wAzgE8AgYHIuXxMozGn3C+C3EdHcYO0K3B8RSwEk3QRsm7dtAdyUG7B1gGdz+Xjgf0hLv38VuLqNa7gTuDEv734icC2wV9mfhNlqzBMfrJo63HOKiHeB54BjgT8B04A9gf6khmJyweSD7SLiOEDAgoLywRGxb8Fh/wTsKalH4alaCeES4NKIGAx8HeiR43oBeEnSXqTG7fdtXMOrEdGc6uhKYGgZH4GZmVVZtWbrTSVlDZ9KapxOJCVafQjYXdIAAEnrS9oWeBzoLWm3XL62pO0Ljvcr4G7gZklrATOBEZI2yUOChxXU3Qj4S359TFFcV5GG934bEe+3FnzR/a6DgIUlX7mZmVVdtRqnacDmwIyIeAl4B5iWh+FGAzdKmk9qrAZGxL+AkcCPJM0jNWSfKTxgRPyUNER4PfASMJY0JPiHXN5sLKkRmwYU31O6A+hJ20N6AKdKWpBjOTXHbGZmddKt0xfl56UujIjh7VYuk9MXmZmVr9T0Rd32OSdJY4BvsGLGnpmZdRHdtnGKiPNJS7d/QNLZrHy/CuDmiDhP0j2kocm1SMOUJ7d1n8rMVub0RVZN3bZxakmemn5eK5v/T0S8oTS3/RZSI/abmgVnZmYfqFluPUkbSJooaV5ec2lUa1kiJB0vaVau+ztJ6+fyw/K+8yRNzWU9JF0tqSlnidgzl4+WdKuke3K2iB+3FV9EvJFfrkV6Xqr73owzM2twtUz8+gXgxYjYMSIGAffQepaIWyNi54jYkTSt+7hcfg6wXy4/KJedDJCfczoCuLbg+aghwChgMDCqtcwPzSRNAl4mrap7Swvbnb7IzKwGatk4NQH7SPqRpOHAlqzIEjGXlOZoi1x3kKRpkppIExqan4GaDlwj6XhSVglIaZKuB4iIRcDzrMgeMSUiXo+Id4A/A1u1FWBE7Ee677QuLWSIcPoiM7PaqNk9p4h4QtJQYH/gh8BkUpaI3Vqofg1wSETMkzQaGJGPcaKkXYEDgLmShpCyTbRmWcHr9ynheiPiHUl3AAfnGM2sBJ74YNVUy3tOfYC3I+IG4AJSSqHWskRsCCzJ2SCOKjhG/4iYGRHnkB643ZKUleKovH1boC8pA0U5sfUsuN+1FqkBXVTxxZqZWYfUcrbeYGCcpOXAu6RnkN4DLpa0UY7lZ8AC4DuklEXPk4YDN8zHGCdpG1JvaQowj9SIXJGHAN8DRucEruXEtgFwh6R1ScOF9wFXdOBazcysA7p1hojO5AwRZmblKzVDRMMt025mZrZaPYQLIGkmaTZeoaMjoqke8ZiZ2aq6TeMkqRdwZERclt+PAE6LiAML60XErq3sP40V97Y+AjwcEYd0XsRmjeOjf5zb4WP8dc8hVYjELOlOw3q9gJMq3TkihjcvfkhamuPWqkVmZmZlaajGSVI/SYskXZXTFE2QtI+k6TkF0S6SxkoaL+l+Sc9IOjXvfj7QX9JcSeNyWU9Jt+RjTlAJU/jycvN7Abd30mWamVk7GnFYbwAp6eoJwCzgSFIWiIOAs0gLEw4kLQW/IfC4pMuBMcCg3PNpHtbbiZRd4kVSdondgQfbOf+hpMwSbxRvkHRCjou+fft25BrNzKwNDdVzyp6NiKaIWE565mlKpPnuTUC/XGdiRCyLiFdIufA2a+VYD0fE4nysuQX7t+UI4MaWNjh9kZlZbTRiz6kw5dDygvfLWRFvqWmJykpfJGkTYBdS78lsteHJDNZoGrHnVKk3WTHbrlKHAXflRLFmZlYn3aZxiohXgel5IsW4dndo2eG0MqRnZma14/RFFXL6IjOz8jl9kZmZdVm1XDLj7pzFodT6/SQ91glx3JafhSr82q+ozlvVPq+ZmZWulosN7l+rc7UlIjwTz1ZrU+7r3ynH3XuvpzvluLZ6qlrPSdIZzdkaJF0o6b78em9JN0h6TtKmuUe0UNKVkhZIulfSernuUEnzJM0ATi449vaSHs69nPmStinIJnFtLrtF0voFx3lA0hxJkwoWEuwv6Z5cPk3SwFz+cUkzJM2S9INqfSZmZlaZag7rTQWG59fDSKmD1iZld5hWVHcb4OcRsT3wGvClXH41cGoLS7efCFyUsz8MAxbn8k8Av4yIHYA3gJPyOS8BRkbEUGA8cF6u/0vgP3L5acBlufwi4PKI2Bn4a2sXKOkESbMlzV66dGn7n4iZmVWkmo3THGBozk23jJQ8dRipwSpunJ6NiLkF+/XLq+H2iogHcvn1BfVnAGdJOhPYKiL+mctfiIjp+fUNpIbwE8AgYLKkucC3gS0k9QQ+A9ycy38BbJ733Z0VU8gLz7sSZ4gwM6uNqt1zioh3JT0HHAv8CZhPyn/XH1hYVL04c8N6pKXXW5zXHhG/zuswHQBMkvQ14JkW6kc+zoLi3pekDwGvNefea+k0bV6gmZnVTLUnREwlDZd9lZQL76fAnIiI9hKCR8Rrkl6XtEdEPAgc1bxN0tbAMxFxcX69A6lx6itpt4iYQcqJ9yDwONC7uTwP820bEQskPSvpsIi4OWco3yEi5pGSwh5O6n0dhVk35okL1hVUeyr5NNJQ2YyIeAl4h1WH9NpyLPDzPCHinwXlo4DH8nDcQOC6XL4QOEbSfGBj0n2jfwEjgR9JmkdK+PqZXP8o4LhcvgA4OJf/X+BkSbOAjcq5YDMzq74umyFCUj9SHrxB9Ti/M0SYmZXPGSLMzKzLasQlM0oSEc+RZuWZmVk30216Tp2V7sjMzGqvy/aczKw0Y8eO7VbnsdVDt+k5ZWsWp0WSdL+kYQA5fdJz+fVoSbdLujNPMT9F0rckPSrpIUkb1/VKzMxWY92tcWotLVJrBgFHkpZmPw94OyJ2ImWk+EpxZacvMjOrje7WOK2SFqmd+n+MiDcjYinwOnBnLm9qaV+nLzIzq43u1jgVp0VaC3iPFdfZo436ywveL8f348zM6mZ1+AX8HDAUeJiUOcJsteKJCtYVdbeeU0suAL4h6U/ApvUOxszM2tdl0xfVm9MXmZmVz+mLzMysy3LjZGZmDWd1mBBhtlpZPKacVWqqZ4vzh9flvNY9dZuek6Rekk4qeD9C0l1l7P8rSfMkzZd0S17W3czM6qDbNE5AL+Ckdmu17j8jYseI2AH4X+CU6oRlZmblaqjGKWcWXyTpKkmPSZogaR9J0yU9KWkXSWMljc85856RdGre/Xygv6S5ksblsp65F7QoH6vVteIj4o0cg4D1gFWmMTp9kZlZbTRU45QNAC4CdiAtyX4ksAdwGnBWrjMQ2I+UE++7ktYGxgBPR8SQiDg919sJ+CawHbA1sHtbJ5Z0NfDXfPxLirc7fZGZWW00YuP5njYXAAASYklEQVT0bEQ0RcRyYAEwJdLDWIX57iZGxLKIeAV4GdislWM9HBGL87Hm0k6uvYg4FugDLARGdfhKzMysIo04W6+UfHct5dBr71ht1ftARLwv6SbgdODqUgI2aySeNWfdQSP2nCr1JrBhJTsqGdD8GvgisKiKsZmZWRkasedUkYh4NU+ceAz4PTCxjN0FXCvpQ/n1POAbnRCmmZmVwLn1KuTcemZm5XNuPTMz67K6zbBeqSTdBny8qPjMiJhUj3jM2vKTUQfWO4SS/ddNJSdkMWvXatc4RcShxWWS1pc0EehPmtV3Z0SMqXlwZmYGeFiv0AURMZD04O7ukv6t3gGZma2uatY4SdpA0sScXPUxSaMkDZX0gKQ5kiZJ2jzXPV7SrFz3d5LWz+WH5X3nSZqay3pIulpSk6RHJe2Zy0dLulXSPTn10Y9biy0i3o6IP+bX/wIeAbZo4RqcvsjMrAZq2XP6AvBiTq46CLiHlCJoZEQMBcYD5+W6t0bEzhGxIylbw3G5/Bxgv1x+UC47GSAiBgNHkKaE98jbhpAyPQwGRknasr0gJfUiPec0pXib0xeZmdVGLe85NQEXSPoRcBfwd2AQMDnnY10TWJLrDpJ0LinTeE+gebLCdOAaSb8Fbs1le5Dz4EXEIknPA9vmbVMi4nUASX8GtgJeaC1ASWsBNwIXR8QzHb5iMzOrSM0ap4h4QtJQYH/gh8BkYEFE7NZC9WuAQyJinqTRwIh8jBMl7QocAMyVNIT00Gxryk1f9EvgyYj4WftXZNb5PAPOVle1vOfUB3g7Im4ALgB2BXpL2i1vX1vS9rn6hsCSnG38qIJj9I+ImRFxDvAKsCUwtbmOpG2BvsDjFcR3LrARKYu5mZnVUS2H9QYD4yQtB94lpQd6D7hY0kY5lp+RMpF/B5gJPE8aDmzOmTdO0jak3tIUUpqhRcAVkpry8UZHxLI2lm5ahaQtgLPzsR7J+14aEVd16IrNzKwiTl9UIacvMjMrn9MXmZlZl7XaZYiQNBNYt6j46Ihoqkc8ZsV+fuJ99Q6hIidfsVe9Q7BupNs0Tvn5pCMj4rL8fgRwWkSslJwsInZtZf8JwDDS/bCHga9HxLudGrSZmbWoOw3r9QJO6sD+E4CBpIkb6wFfq0ZQZmZWvoZqnCT1k7RI0lU5TdEESfvkRQSflLSLpLGSxku6X9Izkk7Nu58P9Jc0V9K4XNZT0i35mBPUxhS+iLg7MlLPyemLzMzqpKEap2wAcBGwA6kncyQpC8RpwFm5zkBgP2AX4Lv5eagxwNMRMSQiTs/1diI9t7QdsDWwe3snz8c6mpReaSVOX2RmVhuN2Dg9GxFNEbGc9MzTlNybaQL65ToTI2JZRLwCvAxs1sqxHo6IxflYcwv2b8tlwNSImNaRizAzs8o14oSIwpRDywveL2dFvKWmJSorfZGk7wK9ga+XGqxZtXnWm1ljNk6VepMVmSTKJulrpKHCvXNPy8zM6qQRh/UqEhGvAtPzRIpx7e6wqitIw4Mz8qSKc6oboZmZlcrpiyrk9EVmZuVz+iIzM+uyGu6ek6R+wF15tdxS6l+T698i6SrgpxHx56I6o4FhEXGKpNuAjxcd5syImIRZHSwc+Ml6h1AVn1y0sN4hWDfScI1TR0REu1kdIuLQWsRiZmaVa9RhvTUlXSlpgaR7Ja0naYikhyTNl3SbpA8X75SzRgzLr4+V9ISkByh4+FbSFyXNlPSopD9I2kzSGjkDRe9cZw1JT0natGZXbGZmH2jUxmkb4OcRsT3wGvAl4DrS8NsOpAdyv9vazpI2B75HapQ+T8oQ0exB4NMRsRPwG+CMPHX8BlasursPMC8/5Ft4XKcvMjOrgUZtnJ6NiLn59RygP9ArIh7IZdcCn21j/12B+yNiaUT8C7ipYNsWwKS8cu7pQPPS8OOBr+TXXwWuLj6o0xeZmdVGo95zKs7s0KuCY7Q2R/4S0qSJO/KyGmMBIuIFSS9J2ovUuB3Vyv5mVeWJBGaratSeU7HXgb9LGp7fHw080Eb9mcAISZvkRK6HFWzbCPhLfn1M0X5XkYb3fhsR73c8bDMzq0Sj9pxacgxwhaT1gWeAY1urGBFLJI0FZgBLgEeANfPmscDNkv4CPMTK08rvIA3nrTKkZ2ZmteMMEQXyTL8LI2J4e3WdIcLMrHylZojoSj2nTiVpDPANfK/JzKzuuso9p04XEedHxFYR8WC9YzEzW911m56TpF7AkRFxWX4/AjgtIg4scf9TSKvm9gd6Fz/jZFZs8LWD6x1CQ2k6pqneIVg30p16Tr2Akzqw/3TSw7fPVyccMzOrVEM1TpL6SVok6aq8LtMESftImp7TC+0iaayk8TlV0TOSTs27nw/0z2sxNa/n1FPSLfmYEySptXNHxKMR8VxnX6OZmbWvEYf1BpCeSzoBmAUcCewBHAScBcwFBgJ7kla+fVzS5cAYYFBEDIEPhvV2ImWAeJHUM9qdlL6oIpJOyHHRt2/fSg9jZmbtaKieU/ZsRDTlfHcLgCmR5rs3Af1ynYkRsSzfF3qZtIJtSx6OiMX5WHML9q+I0xeZmdVGI/acClMXLS94v5wV8RanN2rtOkqtZ1Y2TwAw6zyN2HOq1JukYT4zM+viuk3jFBGvAtPzRIpx7e5QRNKpkhaTspbPz6vqmplZHTh9UYWcvsjMrHylpi/qNj0nMzPrPla7CQKSbmPlTOSQVtidVI94zMxsVatd4xQRh7ZULuk80kq4H46InrWNajUxdqN6R2Cdaezr9Y7AuhEP661wJ7BLvYMwM7MaNk6SNpA0UdK8PKNulKShkh6QNEfSJEmb57rHS5qV6/4uLzCIpMPyvvMkTc1lPSRdLalJ0qOS9szloyXdKumenProx23FFxEPRcSSdq7hBEmzJc1eunRpdT4YMzNbRS17Tl8AXoyIHSNiEHAPcAkwMiKGAuOB83LdWyNi54jYEVgIHJfLzwH2y+UH5bKTASJiMHAEcK2kHnnbEGAUMBgYJWnLjlyAM0SYmdVGLe85NQEXSPoRcBfwd2AQMDnnY12TtKQ6wCBJ55IyjfcEmicrTAeukfRb4NZctgepkSMiFkl6Htg2b5sSEa8DSPozsBXwQqddoZmZVUXNGqeIeELSUGB/4IfAZGBBROzWQvVrgEMiYp6k0cCIfIwTJe0KHADMlTQEaDXTOE5f1Fh8w9zMSlTLe059gLcj4gbgAmBXoLek3fL2tSVtn6tvCCyRtDYFy6ZL6h8RMyPiHOAVYEtganMdSdsCfYHHa3RZZmbWCWrZkxgMjJO0HHgX+AbwHnCxpI1yLD8jZSL/DjCTtPBfEyty5o2TtA2ptzQFmAcsAq6Q1JSPNzoilrWxdFOL8oSJI4H1cxqjqyJibOWXa2ZmlXL6ogo5fZGZWfmcvsjMzLqs1W6CgKSZwLpFxUdHhBfnMTNrEDVvnPLsu3sj4sX8/jlgWF7VtprnuZt0DwngyIi4DCAidm2h7laS5pCms68NXBIRV1QznmL9xkzszMOb1dxz5x9Q7xCsG6nHsN5ooE81DiSp1cY1IvaPiNdIz0qd1M6hlgCfiYghpFmEY/LsQjMzq4N2GydJZ0g6Nb++UNJ9+fXekm6QtK+kGZIekXSzpJ55+zk5BdFjkn6pZCQwDJggaa6k9fJp/iPv3yRpYN5/A0nj8zEelXRwLh+dz3MncK+kzSVNzcd7TNLwXO85SZsC5wP98/YWFyGMiH9FRPMzUeu29rk4fZGZWW2U0nOaCgzPr4cBPfPzR3uQpnl/G9gnIj4FzAa+letemlMQDQLWAw6MiFtynaMiYkhE/DPXfSXvfzlwWi47G7gvInYG9iRNI98gb9sNOCYi9iIN3U3KvZ4dgblF8Y8Bns7nO721i5S0paT5pAwSP2oedizk9EVmZrVRSuM0BxgqaUNSxoUZpEZqOPBPYDvS8uhzgWNIKYIA9pQ0Mz9/tBew/SpHXqE5FdEcoF9+vS9peG0ucD/Qg/SALcDkiPhbfj0LOFbSWGBwRLxZwjWtIiJeiIgdgAHAMZI2q+Q4ZmbWce1OiIiId/OkhWOBPwHzST2Z/sCzpIbiiMJ9cuLVy0gTHV7IDUcPWtc8pFaYYkjAlyJipWwPOX3RPwrimyrps6SURtdLGhcR17V3Xa2JiBclLSA1vrdUepz2+OaxmVnrSp0QMZU03DYVmAacSBo+ewjYXdIAAEnr5xRCzQ3RK/ke1MiCY73JiowPbZlEuhelfOydWqokaSvg5Yi4EvgV8KmiKu2eT9IWzfe/JH0Y2B2nQDIzq5tSG6dpwObAjIh4CXgHmBYRS0mz727M92seAgbmWXJXku5J3U4aemt2DSndUOGEiJb8gDSte76kx/L7lowgJYF9FPgScFHhxoh4lTTs+FhrEyKATwIzJc0DHgAu8HNPZmb14/RFFZK0lJT7r9CmpIS0jahRY3Nc5WnUuKBxY3Nc5ensuLaKiHZnlLlxqiJJs0vJGVUPjRqb4ypPo8YFjRub4ypPo8S1WqUvkjQYuL6oeFlLWSPMzKx+VqvGKd9HGlLvOMzMrG3OSl5dv6x3AG1o1NgcV3kaNS5o3NgcV3kaIi7fczIzs4bjnpOZmTUcN05mZtZw3Dh1gKTDJC2QtFxSq1MvJX1B0uOSnpI0pkaxbSxpsqQn878fbqXej/M1LJR0cXNGjgaIq6+ke3Ncf5bUrxHiynU/JOkvki7tzJhKjUvSEKWVARZImi9pVCfG0+bPsqR1Jd2Ut8/s7O9bmbF9K/8szZc0JWeXqXtcBfVGSoq2fpfUOi5J/yd/Zgsk/boWcX0gIvxV4Rcps8QnSIlph7VSZ03gaWBrYB1gHrBdDWL7MTAmvx5DyrReXOczwPQc45qkpL4j6h1X3nY/8Pn8uiewfiPElbdfBPyalHm/Eb6P2wLb5Nd9SOuT9eqEWNr9WSatnXZFfn04cFNnf0ZlxLZn888R8I1axFbq/39SirWppCw7Lf4uqcPntQ3wKPDh/P4jtfheNn+559QBEbEwihLTtmAX4KmIeCYi/gX8Bji486PjYODa/Ppa4JAW6gQpD+I6pHWs1gZeqndckrYD1oqIyQAR8VZEvF3vuHJsQ4HNgHs7OZ6S44qIJyLiyfz6ReBloDPWdCnlZ7kw3luAvTu7N15qbBHxx4Kfo4eALRohruwHpD9E3qlBTKXGdTzw84j4O0BEvFyj2AAP69XCx0hrRDVbnMs622YRsQQg//uR4goRMQP4I+kv7SWkdbEW1jsuUk/gNUm3Ki00OU7SmvWOS9IawE+AVtcFq0dchSTtQvpj4+lOiKWUn+UP6kTEe8DrwCadEEslsRU6Dvh9p0aUtBtXTmq9ZUTcVYN4So6L9P9wW0nTJT0k6Qs1i47V7CHcSkj6A/DRFjadHRH/U8ohWiiryvz9tmIrcf8BpKHJ5r8gJ0v6bERMrWdcpJ/L4cBOwP8CN5ESDP+qznGdBNwdaRmYjoRS7biaj7M5KQPKMRGxvBqxFZ+ihbLin+VO+3lvR8nnlfRl0pp0n+vUiPLpWij7IK78B8+FpJ/vWirl81qLNLQ3gvQ7YpqkQZESe3c6N07tiIh9OniIxcCWBe+3AFZZZbcSbcUm6SVJm0fEkvxLq6Uu+aHAQxHxVt7n98CnSWPf9YxrMfBoRDyT97k9x9WhxqkKce0GDJd0Euk+2DqS3oqIDk1yqUJcSPoQMBH4dkQ81JF42lDKz3JzncWS1gI2Av5G5yvp/5mkfUiN/uciYlnx9jrEtSEwCLg//8HzUeAOSQdFxOw6xtVc56GIeBd4VtLjpMZqFjXgYb3ONwvYRtLHJa1Dukl8Rw3OewdpZWLyvy318v4X+JyktSStTfpLsrOH9UqJaxbwYUnN9032Av5c77gi4qiI6BsR/Ujrm13X0YapGnHln6vbcjw3d2IspfwsF8Y7Ergv8t30TtZubHn47BfAQTW8f9JmXBHxekRsGhH98s/VQzm+zmyY2o0ru500iQRJm5KG+Z7p5LhWqOXsi+72Rep5LCat5PsS6Z4NpBlTdxfU2x94gnQf4OwaxbYJMAV4Mv+7cS4fBlyVX69J+s+6kPTL/6eNEFd+/3nSqstNpDXA1mmEuArqj6Y2s/VK+T5+GXiXtABo89eQTopnlZ9l4PukX6iQJtjcDDwFPAxs3dmfURmx/SH/P23+jO5ohLiK6t5PDWbrlfh5Cfhp/t3QBBxeq+9lRDh9kZmZNR4P65mZWcNx42RmZg3HjZOZmTUcN05mZtZw3DiZmVnDceNkZmYNx42TmZk1nP8PJo78ZJVvqIQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21218ef2cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE on Training set : 0.39906047938202377\n",
      "RMSE on Test set : 0.7459620508409301\n",
      "r2_score on Training set : 0.8403132358113108\n",
      "r2_score on Test set : 0.6684602242935899\n"
     ]
    }
   ],
   "source": [
    "#### Lasso／L1正则\n",
    "##\n",
    "#   eps（最小alpha除最大alpha的值，默认是0.001，一共100个alpha。也可以像下面一样自己设置alphas）\n",
    "#   max_iter（最大迭代次数）,tol（迭代停止条件，就是那个阈值）\n",
    "#   verbose（是否打出中间结果）\n",
    "#   positive（如果为True会限制系数只能为正数）\n",
    "# class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True, \n",
    "#                                    normalize=False, precompute=’auto’, max_iter=1000, \n",
    "#                                    tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,\n",
    "#                                    positive=False, random_state=None, selection=’cyclic’)\n",
    "\n",
    "lasso = LassoCV()\n",
    "lasso.fit(x_train, y_train)\n",
    "alpha = lasso.alpha_\n",
    "print(\"Best alpha :\", alpha)\n",
    "\n",
    "#画alphas和均方差mse的关系图\n",
    "mses = np.mean(lasso.mse_path_, axis = 1)\n",
    "plt.plot(np.log10(lasso.alphas_), mses) \n",
    "#plt.plot(np.log10(lasso.alphas_)*np.ones(3), [0.3, 0.4, 1.0])\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()    \n",
    "\n",
    "# Plot important coefficients\n",
    "coefs = pd.Series(lasso.coef_, index = x_train.columns)\n",
    "print(\"Lasso picked \" + str(sum(coefs != 0)) + \" features and eliminated the other \" +  \\\n",
    "      str(sum(coefs == 0)) + \" features\")\n",
    "imp_coefs = pd.concat([coefs.sort_values().head(10),\n",
    "                     coefs.sort_values().tail(10)])\n",
    "imp_coefs.plot(kind = \"barh\")\n",
    "plt.title(\"Coefficients in the Lasso Model\")\n",
    "plt.show()\n",
    "\n",
    "y_train_pred = lasso.predict(x_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "\n",
    "\n",
    "y_test_pred = lasso.predict(x_test)\n",
    "y_test_pred += mean_diff\n",
    "\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Test set :\", rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\", r2_score_train)\n",
    "print(\"r2_score on Test set :\", r2_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3  Linear Regression without regularization**\n",
    "最小二乘线性回归\n",
    "最小二乘没有超参数需要调优，直接用全体训练数据训练模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE on Training set : 0.3983895892976355\n",
      "RMSE on Test set : 0.7481497155037375\n",
      "r2_score on Training set : 0.8408497069390944\n",
      "r2_score on Test set : 0.6665127756735301\n"
     ]
    }
   ],
   "source": [
    "# Linear Regression\n",
    "# 线性回归\n",
    "#   fit_intercept(是否拟合截距项),normalize(是否进行归一化，注意标准化还得用StandardScaler，前面做过了这里选false)\n",
    "#   copy_X=True(复制了一份X进行操作，选false会改变原始数据)，n_jobs=1(用多少线程，这里1是单线程，-1是利用所有资源)\n",
    "#class sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)\n",
    "\n",
    "# 1. 生成学习器实例\n",
    "lr = LinearRegression()\n",
    "\n",
    "#2. 在训练集上训练学习器\n",
    "lr.fit(x_train, y_train)\n",
    "\n",
    "#3.训练上测试，得到训练误差，实际任务中这一步不需要\n",
    "# Look at predictions on training and validation set\n",
    "y_train_pred = lr.predict(x_train)\n",
    "y_test_pred = lr.predict(x_test)\n",
    "y_test_pred += mean_diff\n",
    "\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\", rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\", r2_score_train)\n",
    "print(\"r2_score on Test set :\", r2_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 对测试集进行测试，生成提交文件"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据上面实验根据r2_score标准，岭回归和Lasso在测试集上得分分别是0.6703075644943807和0.6684602242935899，岭回归效果好一点，所以用岭回归模型对测试集测试并生成提交文件"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y曾经作过标准化，别忘还原"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test_pred = ridge.predict(x_test)\n",
    "y_test_pred += mean_diff\n",
    "y_test_pred = y_test_pred * y_std +  y_mean\n",
    "\n",
    "df = pd.DataFrame({\"instant\":testID, 'cnt':y_test_pred})\n",
    "#df.reindex(columns=['instant'])\n",
    "#y = pd.Series(data = y_test_pred, name = 'cnt')\n",
    "#df = pd.concat([testID, y], axis = 1, ignore_index=True)\n",
    "df.to_csv('2012_pred.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 366 entries, 365 to 730\n",
      "Data columns (total 2 columns):\n",
      "cnt        366 non-null float64\n",
      "instant    366 non-null int64\n",
      "dtypes: float64(1), int64(1)\n",
      "memory usage: 18.6 KB\n"
     ]
    }
   ],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cnt</th>\n",
       "      <th>instant</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>365</th>\n",
       "      <td>3979.736568</td>\n",
       "      <td>366</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>366</th>\n",
       "      <td>3495.544351</td>\n",
       "      <td>367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>367</th>\n",
       "      <td>3325.810725</td>\n",
       "      <td>368</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>368</th>\n",
       "      <td>3281.387674</td>\n",
       "      <td>369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369</th>\n",
       "      <td>4051.110923</td>\n",
       "      <td>370</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             cnt  instant\n",
       "365  3979.736568      366\n",
       "366  3495.544351      367\n",
       "367  3325.810725      368\n",
       "368  3281.387674      369\n",
       "369  4051.110923      370"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
